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BOTTOM  REFLECTION  OF  UNDERWATER  EXPLOSION  SHOCK  WAVES,  COMPUTER 
PROGRAM 

This  report  is  part  of  a  continuing  study  of  the  interaction  of 
the  underwater  explosion  shock  wave  with  the  ocean  bottom.  The 
computer  program  described  in  this  paper  calculates  the  bottom 
reflection  and  generates  plots  of  the  pressure  history.  Tho  cal¬ 
culations  of  this  program  are  being  used  in  the  bottom  reflection 
study  to  assess  the  potential  danger  to  ships  delivering  nuclear 
underwater  weapons  posed  by  various  bottom  materials. 
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BOTTOM  REFLECTION  OF  UNDERWATER 
EXPLOSION  SHOCK  WAVES,  COMPUTER  PROGRAM 
1.  INTRODUCTION 

The  bottom  reflection  of  the  underwater  explosion  shock  wave 
is  of  interest  to  the  Navy  because  of  the  danger  it  presents  for 
self-damage  to  ships  delivering  nuclear  ASW  weapons.  The  theory 
presently  being  used  to  describe  the  reflection  is  a  linear  spherical 
wave  theory  originally  developed  by  L.  Cagniard  (1)  for  the  calculation 
of  the  reflection  at  an  interface  between  two  elastic  solids.  On 
the  basis  of  Cagniard 's  theory,  Rosenbaum  (2)  derived  equations 
which  describe  the  bottom  reflection  of  underwater  explosion  shock 
waves.  Eritt  (3)  has  greatly  extended  and  generalized  Rosenbaum's 
work.  Britt's  report  should  be  consulted  when  using  the  computer 
program  described  here. 

This  report  describes  a  computer  program,  BOTREF,  written  in 
FORTRAN  IV  for  the  NOL  CDC  6400  computer.  The  code  ca.  _ ulates  the 
pressure  history  of  the  bottom  reflection  of  incident  exponential 
pulses  reflected  from  plane,  homogeneous,  elastic  bottoms  using 
the  spherical  wave  theory.  Major  portions  of  this  program  were 
written  by  the  second  author.  The  first  author  later  brought  this 
program  into  its  present  versatile  form  and  used  it  successfully 
in  practical  applications. 
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The  program  has  options  for  calculating  the  spherical  wave 
reflection  in  two  ways:  (1)  using  real  arithmetic  and  equations 
derived  using  contour  integration  (referred  to  as  the  Cagniard- 
Rosenbaum  method)  and  (2)  using  the  "complex  arithmetic  method". 

The  first  method  is  generally  faster,  but  both  usually  take  less 
than  30  seconds  of  central  processor  time  on  the  CDC  6400  for 
calculating  a  complete  pressure  history.  Also  included  in  the  program 
is  an  option  for  calculating  the  bottom  reflection  using  the 
plane  wave  theory  of  Arons  and  Yennie  (4) .  For  both  the  plane  wave 
and  the  spherical  wave,  the  calculations  include  corrections  for  the 
non-linear  changes  of  the  shock  wave  peak  pressure  and  time  constant 
with  the  distance  from  the  charge. 

The  code  generates  a  CALCOMP  plot  tape  of  the  total  pressure 
history  including  the  incident,  bottom  reflected,  and  acoustic 
surface  reflected  waves.  The  print  out,  in  addition  to  the  pressure 
history,  includes  information  such  as  the  incident  angle,  the  plane 
wave  reflection  coefficient  and  phase  shift,  critical  angles,  arrival 
times,  impulses,  and  energy  flux. 

The  output  of  the  bottom  reflection  program  can  be  directly 
transferred  to  the  PTV  Program  (NOLTR  71-65)  which  is  then  used 
as  a  subroutine.  This  program  calculates  the  jjeak  translational 
velocity  (PTV)  of  a  cylindrical  target.  This  velocity  can  be  used 
as  an  index  for  damage. 

The  equations  used  in  the  BOTREF  code  are  described  in  Section  2 
and  references  are  made  as  to  the  location  in  the  program  where 
each  equation  is  used.  In  Section  3  a  detailed  description  is  given 
of  the  program  organization,  inputs,  outputs,  and  other  important 
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symbols.  The  appendices  contain  a  complete  FORTRAN  listing  of  the 
program,  sample  output,  and  a  CALCOMP  plot. 

The  code  contains  many  comment  cards  so  that  roost  of  the  inputs 
and  outputs  and  much  of  the  organization  is  explained  in  the  program 
listing. 

Comments  on  Terminology.  In  the  acoustic  literature  reflectors 
are  called  either  solids  or  fluids,  depending  on  whether  they  have 
a  shear-strength  or  not.  We  prefer  the  terms  non-rigid  or  rigid, 
because  some  solids,  for  instance,  sand,  have  such  a  low  shear 
strength  that  the  theory  for  a  non-rigid  bottom  yields  sufficiently 
accurate  results,  in  spite  of  the  fact  that  the  material  is  a  solid. 

We  hope  that  our  terminology  will  lead  to  leas  misunderstandings 
than  the  conventional  one  or  the  previously  used  term  "liquid  bottom". 

Rigidity  should  be  understood  as  the  resistance  of  a  body  to 
a  change  in  shape  at  constant  volume.  It  is  equivalent  to  shear 
strength  and  is  measured  either  by  the  Poisson  Ratio  or,  as  in  this 
paper,  by  the  propagation  velocity  of  the  shear  wave.  The  shear 
velocity  is  zero  for  a  non-rigid  material.  Compressibility  is  the 
resistance  to  a  change  in  volume  at  constant  shape  and  can  be 
represented  by  the  propagation  velocity  of  a  compression  wave,  i.e., 
the  sound  velocity. 

The  word  rigid  often  has  the  connotation  of  a  material  having 
infinite  rigidity.  We  use  it  in  the  sense  of  a  material  having  a 
finite,  non- vanishing  rigidity. 
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2.  THEORY  USED  TO  CALCULATE  THE  BOTTOM  REFLECTION 

2.1  Theory  of  the  Bottom  Reflection  of  a  Spherical  Wave 

The  theory  used  in  the  computer  program  described  in  this  report 

has  been  derived  by  Rosenbaum  (2).  Britt  (3)  has  reviewed,  explained, 

and  greatly  extended  Rosenbaum's  work.  A  semi-linear  theory  is 

used  which  describes  all  phenomena  of  interest  with  adequate  accuracy 

The  notation  used  in  this  section  is  essentially  that  of  Britt's 

report.  The  following  exceptions  are  to  be  noted.  We  denote  the 

excess  pressure  by  p  instead  of  P.  Britt  and  Rosenbaum  denote  the 

time  by  t;  we  use  t  for  the  time  and  t  for  Rosenbaum's  reduced 

m 

time  (compare  with  Equation (2.2.2)).  The  program  calculates  the 
step  wave  response  nPm  =  i  Pi  which  corresponds  to  one  reflection 
from  the  bottom.  Multiple  reflections  between  the  surface  and  the 
bottom  are  not  included.  (Multiple  reflections  are  of  minor  impor¬ 
tance  to  underwater  explosion  phenomena  that  lead  to  damage  processes 
When  a  strong  pressure  wave  is  reflected  at  the  water  surface,  most 
of  the  wave  energy  is  left  near  the  surface  and  does  not  propagate 
down  into  the  water  because  of  cavitation  and  spray  formation.) 

We  denote  i  Pi  by  Pr#  the  bottom  reflection  slant  range  i Ri 
by  Rr,  the  incident  or  direct  wave  range  by  R^,  and  the  surface 
reflection  range  by  R  .  We  also  drop  the  subscripts  n  and  m  except 

in  t  and  K  (Equation  2.2.18). 
mm 

The  geometry  of  the  bottom  reflection  is  shown  in  Figure  1. 

The  water  depth  is  H.  The  depths  of  the  charge  and  gauge  are  d 

and  d  .  The  horizontal  distance  between  charge  and  gauge  is  r. 
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we  see  that  the  slant  ranges  are  given  by  the  equations 

Ri  =  “  ^  ) 3  +  r2  ~\  (slant  range  of  incident  wave)  (2.1.1) 

r  3  a  1 1/2 

Rr  »  d  +  r  (slant  range  of  wave  reflected  at  bottom) (2.1.2) 


and 


»  _  I"  /  j  *  ^  a  1  l/a  (slant  range  of  the  wave  re- 

Rs  "  [(dg  +  d}  +  r  J  fleeted  at  the  water  surface),  l2*1*31 

where  dr  =  2H  -  d^  -  d  is  the  depth  of  the  "image"  below  the  gauge. 
Further,  we  have 


and 


cos  9  =  d  /R 
r  r 

sin  9  =  r/R  , 


(2.1.4) 

(2.1.5) 


In  the  water  the  sound  velocity  is  denoted  by  ci  ,  and  the  density 
by  pi  .  Similarly,  the  sound  velocity  in  the  bottom  material  is 
Ca 1  the  shear  wave  propagation  velocity  is  c* ,  and  the  density  is 
Pa .  (Britt  denoted  the  sound  velocity  and  density  of  a  rigid  bottom 
by  c3  and  ps  .) 

2.1.1  Critical  Angles.  For  an  incident  angle  9,  which  is  also 
the  reflected  angle,  the  refracted  or  transmitted  ray  into  the 
bottom  makes  an  angle  9a  (see  Figure  1)  given  by  Snell's  law 


sin  0 


a. 

ca 


sin  0S 


(2.1.6) 


The  angle  83  is  that  angle  at  which  the  pressure  wave  enters  the 
bottom.  Similarly,  the  angle  9*  of  the  shear  wave  in  the  bottom  is 
de  fined 

sin  p  «  g.  Sin  q4 .  (2.1.7) 
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When  c,  >  o,  or  c.  >  cj  the  angles  9a  and  9.  become  90*  at 
incident  angles  9cr  and  9crs  defined  by 


sin  « 


Ci  /ca 


sin  V«=  /c*  • 


(2.1.8) 

(2.1.9) 


ecr  is  called  the  critical  angle  of  the  compression  wave,  and  9 

cri 

ie  called  the  critical  angle  of  the  shear  wave.  These  angles  are 

important  for  calculating  and  interpreting  the  bottom  reflection 
pressure  history. 

2.1.2  The  incident  Pulse.  The  computer  program  assumes  an 
exponential  incident  pulse  p±(t)  given  by 


Pi(t)  =.  PpO^i  exp  [-<t  -  R./ci)/g1  for  t  i  Ri/Ct 
Pi*t*  °  for  t  <  R1/ciJ 


(2.1.10) 


where  G  is  the  time  constant  (usually  denoted  by  °)  and  pp  ,  Pf(R., 

ie  the  peak  pressure  of  the  incident  shock  wave,  a  reduced  notation 

is  used  in  the  machine  program  utilising  the  incident  slant  range 

*l  (Equation  2.1.1)  as  the  characteristic  length.  The  reduced  time  is 

*  -  tc,  (It  is  denoted  by  T  in  the  program).  The  reduced  arrival 

time  of  the  front  of  the  direct  wave  is  thus  t  .  1.  The  incident  pulse 
is  then  given  by 


Pi(t)  =  pp(R.)  exp  [-(t  -  1 )/£f  J  for  t  *  1 

n  i  4-  \  _ 


Pi  (t)  o  0 


for  t  <  1, 


where  g  =  ciG/R^  is  the  reduce^  time  constant. 


(2.1.11) 
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For  the  time  constant  G  and  the  peak  pressure  p^  the  relations 
for  the  actual  underwater  explosion  shock  waves  (high  amplitude 
waves)  are  used  which  when  used  together  with  the  wave  equation 
comprise  the  "semi-linear"  theory.  The  shock  wave  parameters  are 
obtained  from  the  similitude  equations 

G  =  CG  W1/3  (W1/3/^)”0  (2.1.12) 

pp  *  Cp(Wl/3/Ri)np  .  (2.1.13) 

where  CG#  Cp,  nQ,  and  np  are  constants  for  a  given  explosive.  V»  is 
the  charge  weight  in  pounds,  or.  with  appropriate  constants,  the 
yield  in  kilotons.  G  and  p^,  are  calculated  in  the  main  program  in 
Cards  B0TR160-167. 

Examples  of  the  constants  are 


Explosive 

S 

np 

CG 

nG 

TNT 

21600 

1.13 

0.052 

-0.23 

HBX-1 

23800 

1.15 

0.049 

-0.29 

Nuclear 

4. 291* 10® 

1.13 

2.242 

-0.22 

(W  =  Yield 
in  kt) 

4.380* 10® 

1.13 

2.274 

-0.22 

The  values  for  nuclear  explosions  of  the  upper  row  are  the  most 
recent  ones.  Those  in  the  lower  row  are  generally  quoted  in  the  liter¬ 
ature.  The  constants  C  and  C_  are  given  in  psi  and  milliseconds. 

P  G 

2.1.3  The  Surface  Reflection.  The  surface  reflection  p_(t) 

calculated  from  the  simple  acoustic  equation  is 

Ps(t)  »  -Pf(R8)  exp[-(t  -  Rs/ci)/Gs]  for  t  *  Rg/ci  ^ 

p  (t)  =  0  for  t  <  R _/ci  , 

s  s 

where  G  =  G(R„) .  In  reduced  notation  this  becomes 
s  s 
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P,(t)  “  -Pp(R,)  exp  [-(t  -  Rs)/ffs]  for  t  i  Rs 


Ps(t) 


for  t  <  R  , 
s' 


(2.1.15) 


where  g,  -  and  Rg  -  R^.  These  equations  are  coded  in 

Cards  BOTR218,  704,  and  882. 

The  surface  reflection  i,  a  tension  wav.  and  its  pressure  is 

to  be  subtracted  from  the  pree.ure  of  the  incident  and  the  bottom 
reflected  wave. 

Equations  <2.1.14  and  15,  iqnore  cavitation  which  in  sea  wat.r 

0OSS  not  let  pressures  drop  substantial!,  below  the  vapor  pressure. 

in  the  machine  proqram  this  is  taken  into  account  by  a.  teat  that 

makes  sure  that  the  total  pressure  does  not  fall  below  zero  absolute 
(Cards  BOTR713  and  884) . 

for  a  very  oblique  incidence  the  acoustic  treatment  of  the  surface 
reflection  breaks  down  and  must  be  replaced  by  the  anomalous  surface 
reflection  described  in  HOITR  70-31  .  The  machine  program  described 
here  does  not  include  this  mod.  of  the  surface  reflection.  This 
problem  will  be  treated  in  another  machine  program  that  describes 
the  shock  wave  propagation  in  shallow  water. 

2-1,4  ~C°nVOlUtl™  integral.  The  theory  of  the  bottom 

reflection  yields  the  reflected  wave  for  an  incident  step  wave 

This  step  wave  response,  denoted  by  Pr(t,,  iB  the  crucial  point  of 

the  analysis  and  win  be  discussed  in  detail  later,  it  has  the 

dimension  of  (length)'1.  The  pressure  history  of  the  bottom  reflected 

weve  for  an  exponential  incident  wave  pr(t,  is  obtained  from  the 
convolution  integral: 
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Pr(t)  =  Pp(Rr)  £Rr(t)  “  1/Gr  J  e*p[-(t  -s)/gJ  Pr(*)drJ 

6 

for  t  i  6 

(2.1.16) 

Pr(t)  =  0  for  t  <  6  . 

This  equation  is  explained  in  Appendix  D  of  Britt's  report.  The 
scale  factor  p-J,  and  the  time  constant  Gr  are  given  by 

P£  -  RrPF(Rr)  (2.1.17) 

Gr  =  G(Rr),  (2.1.18) 

where  is  the  £lant  range  of  the  reflected  wave,  Equation  (2.1.2) . 
The  factor  Rr  of  the  reduced  pressure  scale  factor  p£  stems  from  the 
definition  of  the  reduced  step  wave  response  Pr(t)  which  includes 
Rr  as  a  factor. 

The  reduced  form  of  the  convolution  integral  is  readily  obtained 
by  introduction  of  "t  =  tci/R^  T,  z,  and  Gr  =  ciGr/Ri. 

The  symbol  in  Equation  (2,1.16)  denotes  the  arrival  time  of 
the  bottom  reflection. 

For  subcritical  incidence,  G  <  9cr*  we  have 

«  -  tc  -  Rr/ci  (2.1.19) 

and  the  reduced  form  is 

6  =  ci  6/Ri  =  Rr/Ri  •  (2.1.20) 

In  this  case  6  is  the  arrival  timetc  of  the  peak  of  the  reflected 
wave . 

For  supercritical  incidence,  6  >  ®crf  the  precursor  of  the 
bottom  reflection  arrives  before  t  =  tc,  namely  at 
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6  =  r/cg  +  dr(Cl"3  -  ca“a)l/s  (2.1.21) 

or  in  the  dimensionless  form 

6  =  rCl/caR.  +  dr[l  -  (ci/ca)2]  l/a/RL.  (2.1.22) 

The  convolution  integral  is  calculated  in  the  BOTREF  program 
Cards  B0TR556,  589,  597,  635,  643,  and  673  using  Simpson's  rule  for 
small  intervals  with  three  equally  spaced  points. 

For  an  exponential  incident  pulse  the  integral  need  not  be 
recalculated  from  t  =  6  for  each  time  step  because  of 

exp(t  +  At)  =  exp(t)  exp(  At). 

The  algorithm  used  to  calcinate  the  integral  in  Equation  (2.1.16), 
which  we  call  Fj, is  as  follows: 

FJ(t)  =  exp(-2At/Gr)FI(t  -  2At)  +{[?r(t  -  2At) exp(-At/Gr) 

+  4Pr(t  -  At) J  exp(-At/Gr)  +  Pr(t)f  At/3  .  (2.1.23) 

This  relation  permits  a  convenient  step-by-step  quadrature  of  the 
integral  using  its  value  for  a  time  2 At  earlier.  The  expression 
is  readily  transformed  into  a  reduced  form  by  the  introduction  of 
t.  At,  and  Gr.  Fj  has  the  dimension  time/length. 

For  supercritical  incidence  Pr(t)  has  a  logarithmic  singularity 
at  t  «  t_.  Since  P  (t)  is  a  rapidly  changing  function  of  t  near 
t  ,  a  smaller  time  increment.  At*  «  At/8,  is  used  in  the  code  for 

w 

the  interval  (tc-aAt,  tc+4At) .  The  code  calculates  a  so  that  there 
are  enough  points  in  the  bottom  reflection  pressure-time  history 
(before  and  after  the  time  increment  change)  to  execute  the  impulse 
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and  energy  flux  integrations.  The  usual  range  is  2.1  <  a  <  6.1. 
Because  these  integrations  are  performed  with  Simpson's  rule  on 
equally  spaced  points ,  each  integration  step  is  completed  on  an 
odd-numbered  point. 

Further,  in  the  time  range  t  -  2At*  <  t  <  t  +  2 At*  we 

c  c 

change  the  integration  variable  in  the  convolution  integral  Fj  to 

v3  »  t3  -  z3  for  z  <:  t 

c  c 

u3  »  z3  -  t3  for  z  i  t  . 

c  c 

The  step  wave  response  Pr(t)  behaves  near  the  singularity  like 

pr«t)  -  Cln(  |  tc  -  t|  ). 

The  change  of  variables  v  and  u  transforms  the  last  two  factors  of 
Equation  (2.1.16)  as  follows: 


P  (z)  dz  =  -  7  P_(z)  dv  z  £  t 

£  «  £  V 


“  Pr(z)  du 


z  a  tc. 


Then  we  obtain 


lim  v  _  ,  . 

z-*t_  -  z  pr(z> 
c 


-  c 


lim  v  c 
z-+t. 


(t»  -  zV7* 


ln(t  -z)  =  0 
c 


2  S  t. 


lim  u  ,  .  _ 

z-»t  z  Pr(z) 

c 


lim 

z-*t 


(z3  -  ty  i/a 


ln(z-tc)  =  0 


z  *  tc. 


This  means  the  integrands  vanish  at  the  singularity  of  Pr,  and 
thus  makes  numerical  integration  possible. 
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Equation  (2.1.24)  below  illustrates  the  variable  change. 

V2Zit' 

pr(t)  =  Pf[pr(t)  -  i)  /ex p[-(t  -  z)/GrJ  Pr(z)  dz 

6 


v<tc) 


“  /  exp[“(t  -  ")/Gr  j  Pr(2)  l  dv 

V  (1^,-2  At') 
u(tc+2At^ 

+  J'  exp  f  -  (t  -  /a  1  d  u  \  i 


«(tc) 


exp  [-(t  -  z)/Gr  j  Pr(z)  |  du 


/  expj-(t  -z)/Grl  Pr(«)  dz  |]. 


(2.1.24) 


V2“' 


Up  to  tine  tc-  24f  and  after  time  tc+  2Af  the  integration 
variable  ia  z  and  the  algorithm  of  Equation  (2.1.23)  i.  used  to 
perform  the  quadrature.  Around  the  eingularity  Simpson's  rule  on 

equally  .paced  interval,  of  v  and  u.  instead  of  z  or  time,  i.  need 
for  the  integration. 

Using  the  algorithms  described  below,  Fj(t)  and  pr(t)  are 
evaluated  in  two  steps  before  and  after  the  singularity.  When 
tc-2At’  <  t  s  tc,  the  following  variables  are  used: 


tl 

=  tc  -  2  At 

(2.1.25) 

Vl 

■  v(ti)  =  (t*-t?)l/s 

(2.1.26) 

ta 

=  {>*  -  (*i/4)8}  l/» 

(2.1.27) 

ta 

"  [fcc  “  (  Vi/2)*]  l/# 

(2.1.28) 
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U  »  It2  -  (v  i/4)a  ]  l/a  (2.1.29) 

The  fifth  time  used  here  is  t  .  However,  P  (t  )  does  not  appear 

c  r  c 

in  the  equations  for  Fj  because  the  transformed  integrand  vanishes. 
The  value  of  at  t  *  ta  is  obtained  from 

Fjtta)  =  Fj  (ti  )exp  [-(t3-tx  )/Gr]  \Pr  (tx  )  exp[-(ta -tx  ) /Gr]vi /ti 

+  3Pr(t3)exp[-(t3-ta)/Gr]vi/t8  +  Pr  (t*  )  Vi  /2t3}  Vi /12  . 

(2.1.30) 

This  equation  is  coded  in  reduced  notation  in  Card  BOTR589.  For 
the  next  step  Fj(tc)  is  calculated  using 

Fl(tc)  ■  FI(ta)exp[-(tc-t3)/Gr]+{  Pr(t3  )  exp[-(tc  -t3  ) /Gr]^/2t3 

+  Pr(t*)  exp  [-(tc-t4)/Gr]  v/t4 }  vx/12  .  (2.1.31) 

This  equation  is  cooed  in  reduced  notation  (Card  B0TR597) . 

Similarly,  after  the  singularity  we  define  the  following 
variables: 


ts  **  t  +  2it*  (2.1.32) 

ux  »  u(te)  *  (tf  -  ta)l/a  (2.1.33) 

c 

ta  a  [ta  +  (ux/4)*]l/3  (2.1.34) 

ta  “  [ta  +  (uj/2)a]l/a  (2.1.35) 

U  -  [ta  +  (3-t1/4)aJl/a  (2.1.36) 


Here  tx  is  the  time  of  the  singularity  tc,  but  Pr(tc)  is  not  needed 
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since  the  transformed  integrand  vanishes.  The  value  of  F^(t«) 
is  then  given  by 

Fjtt*)  -  Fj  ( tc)  exp  [-  ( t*  -tc)  /Gr  ]+  {pr(ts)uiexp  t  (t»-t»)/Gr]/ts 

+  Pr(t.)u1/2t.}u1/12  .  (2.1.37) 

This  equation  is  converted  to  reduced  notation  and  coded  in  Card 
B0TR635.  Then  the  last  step  using  the  special  integration  variables 

is 

Fj(te)  -  Fj  ( t* )  exp  [-  (te  -t8 )  /Gr3  +  {pr  (t*  )  u*exp  C-(te -t»  )/Gr3 /2t* 

+  3Pr(t4)uiexp[-(t6-t*)/Gr3/t4  +  Pr(tB)ui/t6^uj/12  . 

(2.1.38) 


This  equation  in  reduced  form  is  coded  in  Card  B0TR643. 


2.1.5  The  Impulse  and  Energy  Flux.  The  impulse  I  and  energy 
flux  Eg.  are  calculated  in  the  main  program  Cards  B0TR717-766. 

These  calculations  are  made  only  if  the  spherical  wave  bottom 
reflection  is  used.  The  impulse  in  psi-sec  is  evaluated  from  the 
equation 


t 


I 


P(t)  dt, 


where  p(t)  «  p^(t)  +  pr(t)  +  pg(t)  is  the  total  pressure  of  the 

incident,  bottom  reflected. and  surface  reflected  waves  and  t  is 

'  o 

the  time  of  the  beginning  of  the  pressure  pulse  p(t) . 
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Ths  energy  flux  E_  in  in-pei  is  found  from  the  equation 
t  r 

Bf  -  {  /  |p|  P  dt  }  /(2.3066  pid)  . 

*o 

where  2.3066  is  a  conversion  factor  necessary  for  to  be  in  units 
in-psi  when  p  is  in  psi,  time  in  seconds.  Pi  in  gm/cnf ,  and  c»  in 
ft/sec. 

Away  from  the  singularity  of  pr(t)  of  t  ■  tc  and  for  sub- 
critical  bottom  reflections  the  integrals  are  determined  using 
Simpson's  rule  on  equally  spaced  points  as  a  function  of  time. 

Near  the  singularity  the  change  of  integration  variables  is  made 
to  ▼  and  u  as  for  the  convolution  integral.  This  change  of 
variables  is  made  in  Cards  B0TR738-755. 

Also  calculated  in  the  same  section  of  the  program  is  the 
"positive  impulse"  which  is  simply  the  impulse  of  the  positive 
part  of  the  total  pressure  p(t) •  If  the  full  output  option  is  used 
(see  the  input  Z5  in  Section  3.1  and  the  sample  outputs  of  Appendix 
B) ,  the  magnitudes  reduced  impulse  I/W*  ,  reduced  positive  impulse , 
and  reduced  energy  flux  are  calculated  in  Cards  B0TR793-797. 

2.2  The  Cagniard-Rosenbaum  Method  for  Calculating  the  Step  Wave 
Response 

In  this  section  the  Cagniard-Rosenbaum  equations  are  listed, 
and  forms  of  these  equations  similar  to  the  FORTRAN  notation  are 
given.  This  method  is  Sister  than  the  complex  arithmetic  method 
which  will  be  discussed  in  Section  2.3,  but  it  has  the  disadvantage 
that  separate  equations  are  required  for  the  precursor  and  the 
main  wave  and  for  each  type  of  bottom  (determined  by  the  ordering 
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of  ci ,  Ca ,  and  C4).  However,  in  the  coding  we  were  able  to  take 
advantage  of  certain  common  factors  and  terms  and  hence  reduce  the 
number  of  statements  that  would  otherwise  be  required. 

2*2*1  ~n“Riqid  Bottom  Precursors,  a  fast  non-rigid  bottom 
(ca  >  ci)  for  which  9  >  9^  has  a  step  wave  response  at  times 
^  t  <  tc  exPressed  by  the  following  equation  (Britt  (2-1.10)): 


X 

Pr(t)  =  felg  ~  M?  f  1-sin  ^ 

r  J  [(l~b?)  u>a+crabaJ  (tu-N) 


(2.2.1) 


where  cu  =  (o  +  M)/2  +  [  (a  -  m)/2]  sin  tt*/2,  (2.2.2) 

b  =  Pl/P2'  =  «  c,-)l/%  M  -  Vos  9  +  (of*-  0# 

N  =  Tmcos  9  -  (c r3  -  T*)l/a«in  9,  sin  ft  =  r/fc.,  and  cos  »  -  d^. 

In  the  program  the  integration  variable  x  a  tt*/2  is  used.  We ^ Iso 
set,  w  -  cm.  Then  after  rearranging.  Equation  (2.2.1)  can  be  put 
into  the  form  which  is  coded 


where 


2  yfi  b  R.  n/2  f  w  dx 
RiPr<t>  -  -~r— 1  f  — - £  ■ 

r  -n/2  wS  +b2  ( ci  a  aa  -wa ) 


(2.2.3) 


Fx  =  d  -  sin  x){  (c  a+w)  (ci  a-c:M)j/|l+sin  x 
+  4(l-ciT®)l/asin  eA^-ciM)]}*/’, 


(2.2.4) 


-  (1  -  sin  x){[(cos  or  +  w)  P(l)]/  [1*  sin  x+P(2)  ]}*/», 
with  cos  a  ■  cjc  ■  [1  -  (ci/c*)M  1//# , 

P(l)  a  COS  Of  -  ClM, 

P(2)  -  4(1  -  d*^*)*/*  sin  9/p(l) . 

The  variables  cos  a,  P(l) ,  and  P(2)  are  calculated  in  Cards  BOTR238, 
STPA022,  and  23. 
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Th«  integrand  above  ie  evaluated  in  FUNCTION  ONE.  The  variable 
Fx  coded  in  Card  ONE023,  and  the  value  of  the  integrand  is  ONE 
in  Card  ONE055.  The  factor  outside  the  integral  is  calculated  in 
Card  STPA025.  The  integration  for  this  and  all  other  precursors 
is  controlled  by  SUBROUTINE  STPWA  which  uses  the  Gaussian  quadrature 
of  FUNCTION  FGI  to  evaluate  the  integral.  The  value  of  R^Pr(tJ/ 
called  STFW  in  Card  STPA027 ,  is  returned  to  the  main  program  BOTREF 
where  the  convolution  integral  is  executed. 


2.2.2  Rigid  Bottom  Precursor ,  Ca se  ca  >  ci  >  C4 .  The  precursor 
integrands  for  a  rigid  bottom  are  also  evaluated  in  FUNCTION  ONE. 

For  the  case  ca  >  ci  >  c*  (slow  shear)  the  following  equation 


(Britt  (4-1.6))  is  used 

1 


p  (t)  «  f  Ssi^l 

4R_C44  J  (ui-N)1 


L-sin  rci 


(uj-N)1/a  [Aa  +  (B  +  C)a] 


(2.2.5) 


where 

A  =  uj(cTa/2  -  cTa  +  iaa)a  (2.2.6) 

B  -  u)(cTa-wa)  (o8-u)a)1/3|tt)3+cra-cra|  /  (2.2.7) 

C  -  bcT4  (oa-u)a)l/a/4  •  (2.2.8) 

In  a  manner  similar  to  the  non-rigid  bottom,  Equation  (2.2.5) 


can  be  rearranged  to  obtain  the  program  form 


2^2  b  Ri  "/2  FXAX  dx 


(2.2.9) 


where  w  ■  cjuj,  cos  or  «  cia,  and 
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Ax  -  4cj c*4A  -  w[1-2(c4/cx)8  (1-w8) ]8  (2.2.10) 

Bx  ■  4cj c**B  ■  4w(c*/cj  )*  (1-w8)  [  (cosacr 

-  w8)  |(c*/c1)8(w8-l)  +l\V/a,  (2.2.11) 

Cx  -  4CiC44C  -  b(co*a  -  w*)l/8.  (2.2.12) 

As  In  the  previous  case  ?x  and  the  factor  outside  the  Integral  are 
calculated  in  Cards  ONEQ23,  and  STPA025.  The  variables  Ax,  Bx< 
and  Cx  are  coded  in  Cards  ONE043,  044,  and  050.  The  value  of  the 
integrand  is  ONE  in  Card  ONE051,  and  as  before  SUBROUTINE  STPWA 
controls  the  integration. 

2.2.3  Rigid  bottom  Precursor, Case,  ca  >  c,  >  ci .  The  precursor 
for  Ca  >  C4>  ci  (fast  shear)  is  based  on  the  following  equation 
(Britt  (4-2.8)) 

1 

p  (t)  =  f  ?t/2).  at  . 

r  4R  c*4  J  (w-N)1/3  [  A3  +  (B  +  C)3] 

♦l 

H  <2’2’13> 

.  b(o-M)  f  (q+w) 1/3  (A-B)  (l-sin  tt»/2)  dt 

4R  C44  J  (*-N)l/3  [(A  -  B)3  +  C3  ] 
r  -1 

where  ti  =  ^  arcsin  ~  ,.q  "  M  j  . 

(In  Britt's  paper  the  magnitude  B  in  the  second  integral  is  denoted 
by  Ba,  a  precaution  unnecessary  if  the  definition  Equation  (2.2.7) 
is  used.)  Equation  (2.2.13)  can  then  be  written  in  the  form  used 
in  the  program 
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R|Pr(t) 


2y[2  b  R. 


TTR 


r 


TT/2 

/  Fx  Pk  dx' 
-tt/2 


(2.2.14) 


where 

Fk  "  (Ax  “  V/  C(Ax  "  Bx)3  +  Cx8]  for  U)8+ci’a-cTa  <  0  (2.2.15) 
Fk  =  hy  [A*  +  (Bx  +  Cx)a]  for  u)8+c78-cT8  *  0. (2.2.16) 


The  variables  Ax,  Bx,  and  Cx  ara  defined  in  Equations  (2.2.10), 

(2.2.11),  and  (2.2.12).  In  the  first  case  the  integrand  is  coded 
in  Card  OHE047  and  the  second  case  in  Card  0NE051. 

2.2.4  Step  Wave  Response  at  t  =  tc>  At  the  peak  of  the  bottom 

reflection  at  t  =  t  =  R  /ci,the  step  wave  response  P  (t  )  is 

c  i  r  c 

calculated  in  the  main  program  BOTREF.  For  supercritical  incidence, 

9  >  @cr,  Pr  =«  *”  where  the  sign  depends  on  the  phase  shift  c* 
explained  in  Section  2.5.1.  The  treatment  of  this  case  is  discussed 
in  Section  2.1.4.  For  subcritical  incidence,  8  <  0  ,  p  remains  finite 

Cl  4 

and  Pr  ■  K/kr  where  K  is  the  plane  wave  reflection  coefficient  of 
Section  2.5.1. 

2.2.5  Non-Rigid  Bottom  Main  Wave,  Case  cn  >  ci  .  A  fast  non- 
rigid  bottom  (ca  >  ci )  has  a  step  wave  response  at  times  t  >  tc 
given  by  the  equation  (Britt  (2-2.10)) 


pr(t)  *  ± 


1-b 

1+b 


Ik  f - v{°a  ~  [(ou-K  )8+l1  ^  - T ( u)+K  )  8  +l1  }du), 

^r  J  [(l-b^^+aVV  1  m  J  L  m  J  * 


(2.2.17) 
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where 


K  =  t  cos  9 

m  m 


(2.2.18) 


L  =  (Tma  "  c~3) sin3 9 .  (2.2.19) 

The  subscript  m  has  been  kept  to  distinguish  it  from  the  reflection 
coefficient  K. 

In  the  code  the  integration  variable  ie  w  =  ci  to,  and  the  form 

of  the  equation  is  similar  to  that  for  the  precursors 

d  a 


RiPr(t) 


(1-b) R±  2bRi 
( 1+b)  R  +  ~W~ 


f 


Fx  w  dw 


w8+b2 (ci c8  -  wa) 


(2.2.20) 


0 


where  F 

x 


F 


x 


is  now 

r  -  wa  ll/a 

[ci  L+(w-ci  Kj^J 3  J 


_2  2  2  -i/a 

Ci  a  —  w  j 

c8 L  +  (w+ci  Km)  j 


(2.2.21) 


m  f  cos* or  -  w*  i1  _  r  coa*g  -  w*  1 
lP(8)  +(w-P(7)  )*  t>(8)  +(w+P(7)  )*  ^ 

with  cos  a  ■  ci<j  ■  [1  -  (ci/c*)*  l1  . 

The  abbreviations  P (7)  and  P(8)  are  listed  in  Cards  STPB026  and  27. 
The  function  Fx  above  is  calculated  in  Card  ONE032,  and  the  integrand 
is  ONE  in  Card  ONE055.  The  factor  outside  the  integral  is  evaluated 
in  Card  STPB032.  The  first  term  on  the  right  hand  side  of  Equation 
(2.2.20)  is  computed  in  Card  STPB038.  The  integration  for  this  and 
all  other  main  wave  responses  is  controlled  by  SUBROUTINE  STPWB. 

The  value  of  R^Pr(t)  is  calculated  in  Card  STPB047 . 
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2.2.6  Non-Rigid  Bottom  Main  Wave.  Case  ci  >  ca .  The  step  wave 
response  for  a  slow  non-rigid  bottom,  one  with  Ci  >  c? ,  is  expressed  in 
the  equation  (Britt  (2-3.14)) 


Pr(t) 


1  1-b  2j/~2 

R  1^  - 

r 


>i/~2  b  f  via*  -y  )»  /a  <T(u ?  +D)a  +E?  /a+(uf  -F)|^a 

^r  J  ( 1-b8 )  +o* b8  *  (uf+D)8  +  E  * 


where 


and 


5  -  (cT  -  cT*)l/i, 

D  ■  t8  cos  2«  +  cT8  sin8  ft, 
in  > 

E  -  4(sin8ft  cos8 ft)  t^(t^  -  cT3)# 

F  ■  t8  -  cT9  sin?  fl  . 
m 

The  form  used  in  the  program  is_ 

Ci  a 

R  P  ffcl  -  2i.  i=£  _  2%/~2  b  Rj  /*  P  P  dx 

RiPr(t)  R~  I7b  ttr-  ^  FA  FB  dx 

•  ^  A 


(2.2.22) 

(2.2.23) 

(2.2.24) 

(2.2.25) 

(2.2.26) 


(2.2.27) 


where 


Ci  a) 


Fa  =  xtc^o8  -  x8)l/8/[(l-b8)x8  +  ^Ci8^8] 

F  =  cT1  ft33  +  o>3  +  g>X//s  +  (^  -  f)  )  l/a 

B  '  (B*  +  D)8  +  E  ' 


(2.2.28) 

(2.2.29) 


The  integrand  is  evaluated  in  FUNCTION  TWO  Card  TW0017,  Fft  is  coded 
in  Card  TW0013,  and  Ffi  is  coded  in  Cards  TW0014  and  015.  The  terms 
corresponding  to  O,  E,  and  f  are  denoted  by  P(ll),  P(12),  and  P(13) 
and  are  evaluated  in  cards  STPB029-31. 


2.2.7  Rigid  Bottom  Main  Wave,  Case  ca  >  c*  >  ci .  The  rigid 
bottom  main  wave  response  for  the  case  ca  >  C4  >  ci  (fast  shear) 
is  expressed  in  Britt's  equations  (4-4.3),  (4-3.14),  and  (4-3.15) 
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which  are  as  follows: 

Pr(t)  3  R  +  a  (2.2.30) 

r 


+ 


b 

2TTRr<4 


+  2ttr^cJ 


(a*-  (A-B)  |  _ 1 _  _ 1 _  1  d. 

[(A-B)  2+C2  ]  '  [( w-Km) 2  +L]1  's  [( U)+Km)  2  +l]  1  '2' 


.  A  ( o2  -tu8  )l^a 
[  A8  +(B+C)  *] 


{ 


X  X 

[(  --K^)  S+L]l/*  [(ui+K^’+L]*7' 


} 


dw 


where 


=  (Ic4ra 


and 


with 


_  -a/L*  i 

(a8+f)1^8-a(l//a  p 

(2.2.31) 

v>  t 

a8  +  f  ' 

-{»  [<  T 

■  -  k8  )  -  kag3g4j  - 

bga  )  (v  r  04-5  a  8 

-  fc’gjg.l 

-  91  k  [4  ( - k2 ) 

«-  2 

9.  9a  "l 

+  *g»g«^  <*+*>] 

+  .-as-.} 

4g3C4 ) 

• 

Here.  cgt  is  the  propagation  velocity  of  the  Stonley  wave, 
k  »  l/cgfc,  gi  =  (k2  -  cf8)*^8,  93  =  (k2  -  cj2)1^2,  g*  *  (k2  -  cT2)1^2, 
a  a  -  (k2  -  cP  cos8  0),  ard  f  »  4‘rj^g1 8  cos8  9 . 

The  Stonley  wave  propagation  velocity  cgt  is  calculated  in 
SUBROUTINE  STONL.  The  equation  for  c  ^  used  in  the  program  is 
described  in  Section  2.4. 
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The  above  equation  is  coded  in  the  form 
R,  2bR.  Cl/ 

RiPr(t)  “  R  +  RiA  +  - "  /  pxFk  ^  ‘  (2.2.32) 

r  0 

Fx  and  have  been  defined  in  Equations  (2.2.21),  (2.2.15),  and 
(2.2.16). 

The  first  two  terms  of  Equation  (2.2.32)  are  calculated  in 
Cards  STPB056-71  for  all  solid  bottom  main  waves,  and  the  result 
is  stored  in  the  variable  TERMl.  The  integrand  is  determined  in 
FUNCTION  ONE  in  Cards  ONE047  and  051  in  the  same  way  as  for  the 
precursor.  However,  the  function  Fx  and  the  factor  in  front  of  the 
integral  are  here  calculated  in  Cards  ONE032  and  STPB032  as  they 
were  for  a  fast  fluid  bottom  main  wave. 

2.2.8  Rigid  Bottom  Main  Wave,  Case  ca  >  ci  >  C4 .  The  main 
wave  response  for  the  rigid  bottom  case  ca  >  ci  >  c*  (slow  shear) 
is  given  by  the  following  equation  (Britt  (4-3. 13)) 


Pr(t) 


+  A 


+ 


_b_  J-W4SJ Ul*  i  —1 -  ,  A - i - ,  }  d. 

2TTRrC*V  A8  +  (B+C) 8  t[(tt)-Kjn)a+L]  [(a)+Kin)a+L]1/aj 


yfi  b  r  (&+ns)l/BK  f  +  csP-f)  j1  /a 

2ttR_c*4  J  (A+C)a+B8  '  (HP  +D)°  +E  * 

r  0 


where 


A  -  uj  [ c*8 /2  -  cT8  -  uj*  ] 


8 

> 


(2.2.33) 


(2.2.34) 
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B  ■  t»  (cT*  +  ui*)  (o8  +  JSP)1  (ca®  -  tf)l//8  , 


(2.2.35) 


C  -  — (o8  +  (S8)1/*. 

4c*4 

The  above  equation  is  coded  in  the  form 

ci  a  ci  aa 


Vr<*>  -  +  V  +  V*.  **  -(|-f  4/'¥b 

r  o  o 

where  x  ■  ci  w,  Fx  and  are  defined  in  Equations  (2.2.21) 
and  (2.2.16), 


(2.2.36) 

dx|  (2.2.37) 
,  (2.2.15), 


(Sf  +  o8)*/8  §  (x8  4-  cos8flr)1//8Bx 

Ci4[(A-W)a  +  B8]  (Ax  +  c^8  +  fi8 

l  (Rtf  +  to)8  +  Bl1/8  +  (tf  -  F)  V1 
1  (E8  +  D)a  +  E  >  ' 


(2.2.38) 


(2 .2  3  39) 


cos  a  ■  CiC  *  Cl  —  (Ci/Cfc,)*)  1  f%, 

Ax  -  Ci® A  -  x[  (ci/c*)®/2  -  1  -  x8  ]  •  , 

Bx  "  Cl®  B  *  X(1  +  X*)  {[cOS®Of+jfl  [(Cj  /c*  ) ®-l-X*l}  1  ^®  , 

Cx  ■  ci®C  ■  b(ci/d*)4  (cos*a  +  x*)x^®/4  . 

The  first  three  terms  of  Equation  (5*2.37)  are  calculated  using  the 
seme  cards  as  for  the  fast  shear  case.  The  integrand  of  the  second 
integral  is  computed  in  FUNCTION  ONE1.  fa  and  Ffi  are  expressed  in  Cards 
0NE1019,  20,  and  21.  A.  ,  §x,  and  Cx  are  calculated  in  Cards  ONE1015-17. 
The  terms  corresponding  to  D,  E,  and  F  are  denoted  by  P(ll) ,  P(12) ,  and 

P(13)  and  are  evaluated  in  Cards  STPB029-31.  The  value  of  the  integrand 
is  stored  in  the  variable  0NE1  in  Card  ONE 102 3.  The  response 
STPW  •  R^pr (t)  is  then  determined  in  Cards  STPB079  and  80, 
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2.3  The  Complex  Arithmetic  Method  for  Calculating  the  Step  Wave 
Response 

A  second  option  for  calculating  the  step  wave  response  is 
provided  by  the  complex  arithmetic  method.  This  procedure  is  based 
on  the  equation  (Britt  (5-2.12)) 

y a 

P r<t}  “  W  /  He  (K  -  %)}  dy 

yi 


2  (  r  H  <Ma  -  d  t/fe 


(2.3.1) 


where  u 


x  +  iy  and  for  t  <  t  =  R  /ci 

c  r 

x  =  0 


yi  =  cp 


ya  =  [tr  -  dr(cT»Rj  -  ta)  l/a] . 


(2.3.2) 


For  times  t  >  t  these  variables  are 
c 


x  3  Rr*  dr*tS  "  1  ^ 


Yi  -  0 

ya  =  Rj?  tr  o 

The  reflection  coefficient  K  for  a  solid  bottom  is  defined 


(2.3.3) 


<*1  []  (2u8  +c?)2  -  4ua  ey4]  -  baaC^4 

K  =  -  ,  (2.3.4) 

o' i  [  (2ua  tcT3)8  -  4ua  or*  J  +  b'VacT4 

where  (v^  =  (cTs  +  ua)1//a  for  i  =  1,  2,  4. 
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For  a  fluid  bottom  c*  =0,  and  the  equation  for  K  reduces  to 

K  =  (fVi  -  hota)/(a\  +  ba9)  .  (2.3.5) 

Ki  is  the  value  of  K  at  u  =  x  +  iy3 .  The  other  variables  used  above 
are  as  follows: 


Y  *=  [u3r8  +  (t  -  drcvi)8J  1 
on  »  £  ci"2  +  (x  +  iyi )®  3  1  ^2 
o*  =  [cp  +  (x  +  iya)2]  1/2 

f  ( )  =  [_  R2 uj!  -  2drtu)a+  (ta-cr2i?)]  1/2 


(2.3.6) 

(2.3.7) 
(2.3.8) 

+  tt*Rr  -  (2.3.9) 


The  form  of  Equation  (2.3.1)  which  is  coded  is 

/  ca  y?  >  TO 

RiPr(t)  ={A2  +  f  W* [F •<*-*)]  dz^/(~*) 

<5yi  i 

where  z  =  Cjy, 


(2.3.10) 


Aa  =  Im[l<!  log  j[ci  dg  -Ci  t^cosoJ  /[(ci  8  coj  8  -2ci  t^cosQ  (Ci  iu,  ) 

+  (ci®Tm=^  sinafl))l/8  +  Clu4  -  q^cos  e]|J^  (2.3.11) 


and 


Cl  u 


UR 


F  = 


ci  ^  (c,  Y/Rr)  ci^y 


(2.3.12) 


As  in  the  Cagniard-Rosenbaum  method,  the  response  STFW  =  R^Pr(t) 
is  calculated  in  SUBROUTINE  STPWA  for  the  precursor  (t  <  t  )  and 
in  SUBROUTINE  STPWB  for  the  main  wave  (t  >  tc)  using  the  Gaussian 


quadrature  of  FGI  to  evaluate  the  integral.  The  last  factor  in 
Equation  (2.3.10)  is  calculated  in  Cards  STPA039  and  STPB097.  The 
integrand  and  Aa  are  coded  in  FUNCTION  SEVEN,  Cards  SEVN035  and  045. 
The  value  of  A8 (t)  is  obtained  from  STEWA  by  a  call  to  SEVEN  with 
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z  =  ci  ya .  The  function  Ki  =  K(x+iya)  is  evaluated  using  the  same 
equations  as  for  K(u)  in  the  integral,  namely,  RCOE  in  Cards  SEVN022 
and  029.  m  FUNCTION  SEVEN  the  variables  brought  over  by  COMMON 
statements  are  calculated  in  the  main  program,  and  members  of  the  P 
array  are  determined  in  Cards  STPA035-38  for  the  precursor  and  in 
Cards  STPB093-96  for  the  main  wave. 


2.4  The  Stonley  Wave  Propagation  velocity 

The  Stonley  wave  propagation  velocity  c  is  defined  as  the 

St- 

zeroes  u  =  ±  i/cgt  =  denominator  of  the  solid  bottom 

reflection  coefficient  expressed  in  Equation  (2.3.4).  Thus 
ua  »  -c~t  is  the  solution  of  the  equation 


a4 

*  (2ua  +  cJT8)2  - 

4ua  a4J  +  bargcT4  ~  0, 

(2.4.1) 

where 

®i  *  (c T3 

+  ua)l/a. 

ora  =  (C P 

+  u8)l/a, 

and 

*4  *  (cT8 

+  ua)l/a. 

To  obtain  the  form  of  Equation  (2.4.1)  which  is  used  in  the  program, 
first  note  that  the  square  roots  o-i  ,  as ,  and  ar4  are  imaginary  since 
cgt  is  known  to  be  smaller  than  cj  ,  ca  ,  and  c*  .  Next  replace  u8 
by  -c“8t  multiply  through  by  icicaC44Cgt,  and  set  ya  =  c8tto  obtain 


/a 


(Cia-  ya  )  |ca(ya  -  2C48)8  -  4c43  [(c88  -  ya)(c48  -  y8)]1//af 

(2.4.2) 

+  bcj ya 8  ( ca 8  -  ya)1/s  =  0. 
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This  equation  is  solved  for  y2  in  SUBROUTINE  STONL  by  iteration 
using  the  secant  method.  The  variable  y3  is  denoted  by  the  FORTRAN 
symbol  Y2.  Then  cgt,  called  CSTON  in  the  code,  is  the  square  root  cf 


2.5  Theory  of  the  Plane  Wave  Bottom  Reflection 

In  cases  where  the  plane  wave  bottom  reflection  is  adequate 
for  one's  needs  or  when  one  wishes  to  compare  these  results  with 
the  spherical  wave  reflection,  the  plane  wave  option  of  the  BOTREF 
program  can  be  used.  The  reflection  geometry,  incident  and  critical 
angles,  and  the  incident  pulses  are  the  same  as  for  the  spherical 
wave  in  Section  2.1;  and,  unless  otherwise  noted,  the  notation  is 
the  same . 


2.5.1  The  Plane  Wave  Reflection  Coefficient  and  Phase  Shift. 
The  plane  wave  reflection  coefficient  K  and  phase  shift  <t>  for  a 
non-rigid  bottom  are  calculated  from  the  following  equations.  For 


subcritical  angles  of  incidence  K  and  5  are 

K  =  (AT  -  1)/(At  +  1) 


(2.5.1) 


<5  =  0 


where 


A  =  cos  9/  [  b (sin2 0  -  sin2^)1^2  ].  (2.5.2) 


At  the  critical  angle  ncr  these  expressions  reduce  to  K  =  1 
and  <5=0.  At  supercritical  incidence  we  have 


I  K  \  =  1 


(5  =  2  arctan  [  b(sin2  t  -  sin2 )l  /cos  e]. 


(2.5.3) 
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The  above  equations  are  coded  in  the  main  program  Cards  B0TR260-274, 
307.  The  FORTRAN  variables  CR  and  E2  denote  X  and  0.  If  K  is 
complex,  then  CR  =  |  K  |  . 

For  a  rigid  bottom  K  and  4  are  determined  from  the  equations 

below.  At  subcritical  incidence,  0  <  8  <8  ,  we  have 

cr  crs 


+  bt  -  1)/(At  +  bt  +  1) 


0  a,  0 


(2.5.4) 


where 


At  =  cos  e  [l  -  2sin8°/sinaQcr&]a/[b(8inaocr-sinae)1 

(2.5.5) 


B_  =  4cos  0sina9(sina0  -  sin8**)/ 
x  cr  8 

[b  sin4ecrs(sinaecrs  -  sina8)l/a] 


(2.5.6) 


At  the  critical  angle  0  =  ©cr  the  equations  simplify  to  K  =  1  and 

0  =  0.  For  an  incident  angle  in  the  range  8  <  9  <  9  the 

cr  wir  s 

reflection  coefficient  is  complex.  Its  modulus  is 

|K|  =  {[\A  +  (BT  -  f  (BT  +  1)8  (2.5.7) 


and  the  phase  shift  is 


where 


0  »  arctan[(l-BT)  ApA]  +  arctan  [(1+BT)ATA]  ,  (2.5.8) 

Ata  =  cos  9  £  l-2sina  9/sina 9 crsJ3 /  £b(sina°-sinaD^r)1 /^J(2.5.9), 


At  the  critical  angle  of  the  shear  wave  8  the  equations  reduce  to 

cr  s 

IKI  =  1 
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and 

0=2  arctan(l/ATA) . 

For  angles 

of  incidence  e  >  0  we  have 

crs 

|K|  -  1 

and 

«  *  2  arctai[l/(ATA  +  BTA)j 

where 

Bta  =4  cos  0  sin8 9 (sin2 «crg  -  sin8 ft)/ 

[b  sin4ftcrs(sin8ft  -  sin8ftcrg)l/a 

(2.5.10) 


(2.5.11) 

(2.5.12) 


These  equations  for  the  solid  bottom  reflection  coefficient  and 
phase  shift  are  coded  in  Cards  B0TR280-3C7.  As  for  the  fluid  bottom, 
K  and  0  are  denoted  by  CR  and  E2;  and  if  K  is  complex,  CR=  1KI. 

2.5.2  The  Plane  Wave  Bottom  Reflection  Pressure  History. 

The  plane  wave  bottom  reflection  pressure  history  pr(t)  is  calculated 
from  the  following  equations: 
when  ft  s  6 _ j 


Pr  =  0 


Pr  =  P^(Rr)  K  exp[-  (t-tc)  /Gr] 


for  t  <  t  a* 


for  t  ^  t 


Rr/Cj 


(2.5.13) 


when  0  >  3 


cr  ' 


pr  =  pF(RrJ  ’|,exp[-(t-tc)/ar]Ei[(tc-t)/Gr]sin  0  for  6  *  t  <  t 
Pr  =  ±  “  with  the  sign  of  0 

pr  =  PF(V  iKlexpf-ft-t^/Gj^  os  0 


(2.5.14) 
for  t  =  t 


-  ^  Ei  j(t-tc)/Gr  sin  0^ 


for  t  >  tc*  (2.5. 


Note  that  the  plane  wave  theory  has  been  modified  to  use  p  (R  ) 
and  Gr  =  G (Rr)  which  account  for  non-linear  changes  of  the  shock 
wave  peak  pressure  and  time  constant  with  distance.  Also  the  arrival 
times  of  the  main  wave  and  precursor  have  been  changed  to  conform 
to  the  sphericall  wave  situation.  In  the  strict  plane  wave  theory 


15) 
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the  precursor  begins  at  t  =  -  00 ,  and  the  incident  wave  and  the 
reflected  peak  arrive  simultaneously. 

The  functions  Ej  (x)  and  Ei(x)  are  the  exponential  integrals 
defined  for  x  >  0  as 

GO 

Ei  (x)  =  /  dy  (2.5.16) 

X 

oo 

Ei(x)  =>  -  /  exPHd  dy  -  -Ex  (-x) 

j  y 

-X 

X 

=  j  dy  .  (2.5.17) 


The  function  Ei  (x)  is  evaluated  using  the  following  approximate 

formula  (see  for  example  Abramowitz  and  Stegun  (5)) 

0  *  x  <  1 

Ei  (x)  a0  +  ai  x  +  a2xs  +  agx3  +  a4x4  +  agx?  -  log(x) 

(2.5.18) 


a0 

=  -.57721566 

a3  =  .05519968 

ai 

=  .99999193 

a4  =  -.00976004 

a3 

=  -.24991055 

a5  =  .00107857 

x  £  1 

x4  +  ai  x3  f  a2  x2  +  a3  x  +  a4 
x  exp(x)  Ea  (x)  F&  ”  >-■■■<  ■  ■  ■  ,.-■  ■■  — --  ■■■  —  ■ 

x4  +  bix3  +bsx2  +  b3X+b4 


(2 .5.19) 


ai  =  8.5733287 

a2  =  18.059017 


bi  =  9.5733223 

ba  =  25.632956 
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a3  =  8.6347609  i*  =  21.0996531 

a4  =  .26777373  b*  =  3.9584969 

The  function  Ei(x)  is  evaluated  for  x  *  .5  using  the  formula  (Reference 
(5)) 

’  xn 

Ei(x)  »  Y  +  log  (x)  +  2^  (2.5.20) 

where  y  -  .57721566...  is  Euler's  constant.  For  x  >  .5,  Ei(x)  is 
obtained  from 

exp(-x)Ei(x)  =  exp(-x)Ei(l)  +  J  dy,  (2.5.21) 

1 

where  E  i(l)  =  1.8951178.  The  integral  is  then  evaluated  using 
the  Gaussian  quadrature  of  FUNCTION  FGI. 

The  reflected  pressure  pr  =  FBOT  is  calculated  in  the  main 
program  in  Cards  BOTR861-879.  The  exponential  integrals  Ej  and  Ei 
are  calculated  in  the  subprograms  EXE1  and  EXEI,  and  the  integrand 
of  Equation  (2.5.21)  is  coded  in  FUNCTION  EXPO. 
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3.  THE  BOTTOM  REFLECTION  COMPUTER  CODE 

The  Bottca  Reflection  Code  has  been  programmed  in  FORTRAN  IV 
for  use  on  the  CDC  6400  computer  at  NOL.  The  code  is  made  up  of 
a  main  program  called  BOTREf  and  the  following  bottom  reflection 
related  subprograms:  STONL,  STPWA,  STPWB,  ONE,  ONE1,  TWO,  SEVEN, 

EXE1,  EXEI,  EXPO,  FGI,  PLOTl,  and  SCAL.  In  addition,  the  NOL 
general  purpose  plotting  program  CALCMl  must  be  included  for  the 
generation  of  a  tape  to  be  plotted  on  CALCOMP  incremental  plotters. 
For  NOL  users  CALCMl  is  available  on  the  subroutine  library  tape. 

The  control  cards  which  are  included  in  the  program  listing  of 
Appendix  A  contain  the  statements  necessary  for  using  CALCMl  from 
the  library  tape.  For  programmers  outside  of  NOL  information  on 
the  plotting  programs  may  be  obtained  from  the  NOL  Mathematics 
Department  (Code  330) . 

The  basic  organization  of  the  bottom  reflection  code  is  as 
follows.  The  main  program  BOTREF  handles  all  of  the  input  and  output 
and  calculates  the  shock  wave  peak  pressure  and  time  constant  and 
other  time  independent  magnitudes.  It  performs  the  time  incrementa¬ 
tion  of  the  pressure-time  histories  and  calculates  the  convolution 
integral,  impulse,  and  energy  flux  for  the  spherical  wave  bottom 
reflection. 

The  spherical  wave  step  wave  response  Pr(t)  is  obtained  by 

calls  from  BOTREF  to  STPWA  for  the  precursor  and  to  STPWB  for  the 

main  wave.  These  subroutines  in  turn  set  up  the  integration  for 

Pr(t)  using  the  Gaussian  quadrature  in  FGI.  The  various  integrands 

described  in  Sections  2.2  and  2.3  are  calculated  in  subprograms 

ONE,  0NE1,  TWO,  and  SEVEN.  The  Stonley  wave  propagation  velocity 
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c  t  for  rigid  bottoms  is  computed  in  SUBROUTINE  STONL  on  a  call 
from  the  main  program. 

The  plane  wave  bottom  reflection  is  also  calculated  in  the  main 
program.  Calls  to  SUBROUTINES  EXE1  and  EXEI  obtain  the  exponential 
integrals  E*  and  Ei  which  are  used  to  determine  the  bottom  reflection 
in  Equations  (2.5.13),  (2.5.14),  and  (2.5.15). 

SUBROUTINES  PL0T1  and  SCAL  set  up  the  CALCOMP  plots  of  the 
pressure-time  history.  PLOTl  calls  SCAL  to  scale  the  plot,  calls 
CALCMl  for  plotting  the  axes  and  the  pressure-time  curves,  and  then 
calls  SUBROUTINES  SYMBL4  and  NUMBR,  which  are  part  of  the  CALCMl 
program,  to  write  additional  information  on  the  plots. 

The  Bottom  Reflection  Program  also  has  an  option  for  calculating 
the  peak  translational  velocity  (PTV)  induced  in  submerged  or  floating 
targets  by  the  bottom  reflected  pulse.  Either  of  the  spherical  or 
plane  wave  reflection  theories  may  be  used.  The  targets  are  approxi¬ 
mated  by  an  infinitely  long  cylinder  of  a  specified  radius,  and  the 
PTV  Program  described  in  Reference  (6),  is  used  to  calculate  the 
peak  translational  velocity.  This  program  uses  the  additional 
subroutines  PTV,  FV,  FI,  XMAX,  VTAB,  and  PTAB.  The  PTV  is  calculated 
by  calling  SUBROUTINE  PTV  (Cards  BOTR813L  and  813R)  . 

The  cards  in  the  main  program  which  are  necessary  for  PTV 
calculations  are  denoted  by  card  numbers  followed  by  letters  A,  B, 

C,  etc.  If  the  bottom  reflection  program  is  not  to  be  used  for 
PTV  calculations,  these  cards  and  the  subroutines  of  the  PTV  Program 
may  be  omitted. 

In  the  following  paragraphs  the  most  important  FORTRAN  symbols 
of  each  subprogram  are  described,  and  the  locations  in  the  program  are 

given  where  each  symbol  is  calculated. 
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3.1  FORTRAN  Symbols  of  the  Main  Program 
Program  Input 

The  input  data  is  read  in  Statements  3  and  4,  Cards  B0TR041,  42, 
and  89,  and  in  Card  BOTRlOlI  using  the  format  8F10.5.  These  inputs  are 
explained  in  comment  Cards  BOTROll-39,  72-87,  and  10JLB-101G.  The 
inputs  and  their  units  are  as  follows: 

First  Data  Card,  Statement  3 
WCH  charge  weight  W  in  pounds  or  KT 

CWAtER  sound  velocity  ci  of  water  in  ft/sec 

CBOT  sound  velocity  c3  of  the  bottom  material  in  ft/eec 

CSHEAR  a  double  purpose  input  expressing  the  rigidity  of  the  bottom. 

If  CSHEAR  >  .5,  it  is  the  shear  wave  propagation  velocity 
c«  of  the  bottom  in  ft/sec.  If  CSHEAR  *  .5,  it  is  the  dimen¬ 
sionless  Poisson  ratio  from  which  the  shear  velocity  c*  is 
calculated  in  Card  B0TR062.  values  of  c*  «  .5  can  be  neglected. 
RHOWAT  density  Pi  of  water  in  gm/cm3 
RHOBGT  density  p3  of  the  bottom  material  in  gm/cm3 
PRECOE  coefficient  C^  of  the  pressure  similitude  equation  in  psi  . 

PRECOE  depends  on  whether  W  is  in  pounds  or  KT. 

Z5  a  control  parameter.  Z5  greater  than  zero  results  in  a 

shorter  print  out  for  the  spherical  wave  reflection. 

See  Appendix  B  to  compare  the  short  and  long  print  out . 

Second  Data  Card,  Statement  3 

PREEXP  exponent  n^  of  the  pressure  similitude  equation 
THECOE  coefficient  CG  of  the  time  constant  similitude  equation  in 
seconds.  This  variable  also  depends  on  the  units  of  W. 
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THEEXP  exponent  n_  of  the  time  constant  similitude  equation 

STEPS  number  of  points  in  the  pressure-time  history  for  one  time 
constant  G.  STEPS  =  20.  is  usually  sufficient  to  obtain 
a  smooth,  detailed  pressure  history.  In  many  cases, 

STEPS  =  10.  or  5.0  is  adequate. 

DURAT  duration  of  pressure-time  history  in  multiples  of  the  time 

constant  G.  If  negative,  its  absolute  value  is  the  duration 
after  the  arrival  of  the  bottom  reflection  peak  at  t  ®  t  , 

If  positive,  it  is  the  duration  after  the  direct  wave 
arrival.  DURAT  =  -3.0  is  generally  sufficient  for  calculating 
the  significant  parts  of  the  bottom  reflection. 

XI  CALCOMP  plot  scaling  parameter  for  the  Y-axis  in  psi  per 

inch  of  graph.  The  X-axis  is  always  drawn  three  inches 
above  the  bottom  of  the  graph.  The  length  of  the  Y-axis 
is  nine  inches.  Thus  the  maximum  pressure  plotted  is 
6  *  XI,  and  the  minimum  is  -3  *  XI.  Pressures  outside  of 
this  range  are  plotted  at  the  maximum  or  minimum,  whichever 
is  applicable. 

X2  scaling  parameter  for  the  X-axis  in  microseconds  per  inch 

of  graph.  If  X2  0.,  SUBROUTINE  SCAL  calculates  an 
appropriate  value  of  X2. 

SLOPE  slope  of  the  bottom  from  charge  to  gauge  in  degrees.  If 
the  slope  is  not  zero,  the  internal  computing  geometry  is 
changed  in  Cards  BOTR170-183.  SLOPE  must  be  zero  if  the 
geometry  changing  options  of  Z2  and  THOVAL  are  used. 

Third  Data  Card,  Statement  4 

BIGH  water  depth  H  at  the  charge  in  feet.  BIGH  is  also  used  as 
a  control  parameter.  After  completion  of  each  bettom 
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reflection  pressure-history,  the  program  control  returns 
to  Statement  4  to  read  a  new  set  of  data.  If  BIGH  =  0., 
the  program  stops.  If  positive,  computation  continues  with 
the  new  geometry.  If  negative,  program  control  transfers 
to  Statement  3  where  a  new  set  of  charge,  physical  constants, 
etc*,  are  read. 

D  depth  d  of  the  charge  below  the  water  surfac^Kn  feet 

DGAU  depth  d^  of  the  gauge  in  feet  *1 

SMALLR  horizontal  range  r  between  charge  and  gauge  i®feet 

TEOVAL  desired  ratio  between  the  bottom  reflection  incident  angle 

0  and  the  critical  angle  9  <,  The  variables  D  and  DGAU 

cr 

are  changed  in  Cards  S0TR137-142  to  obtain  this  ratio. 

SMALLR  is  not  changed.  If  THOVAL  *  0  the  geometry  is  not 
changed.  See  Appendix  C  for  a  discussion  of  this  option. 

Zl  parameter  which  selects  the  theory.  When  Z1  =  0.  the 

spherical  wave  Cagniard-Rosenbaum  method  is  used.  When 
Zl  =s  1.0,  the  Arons-Yennie  plane  wave  theory  is  used.  And 
for  Zl  =  3.0,  the  complex  arithmetic  method  is  used  to 
calculate  the  spherical  wave  bottom  reflection.  Cards 
B0TR389-4^3  make  the  theory  selection  and  write  out  the 
appropriate  headings. 

Z2  arrival  time  difference  between  the  bottom  reflection  peak 

(at  t  =  t  )  and  the  direct  wave  in  microseconds.  If  Z2  £  0., 
c 

the  geometry  is  not  changed.  When  the  geometry  is  changed, 

D  and  DGAU  are  varied  to  obtain  the  desired  arrival  time 
difference.  SMALLR  is  not  changed,  and  the  change  in  D 
is  the  negative  of  the  change  in  DGAU  so  that  the  incident 
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angle  0  is  also  unchanged.  This  geometry  change  is  performed 
in  Cards  BOTR121-127.  See  Appendix  C  for  a  discussion  of  this 
option. 

Z3  plot  control  parameter.  A  CALCOMP  plot  tape  is  generated 

if  Z3  =  0. 

Fourth  Data  Card  (BOTRlOlI) ,  For  PTV  Calculation 
RADIUS  cylinder  radius  in  feet.  This  is  the  draft  or  cross- 
sectional  radius  of  the  target  vessel.  If  RADIUS  £  0o, 
the  PTV  is  not  calculated. 

APRINT  controls  printing  in  SUBROUTINE  PTV.  If  APRINT  *  0., 

the  translational  velocities  calculated  in  the  iteration 
for  the  PTV  are  printed.  An  example  of  this  printout  is 
given  in  Table  B.l  following  the  pressure-time  history. 

If  APRINT  *  o.f  the  variables  TIME1,  PTVl,  and  PTV2  described 
below  are  printed  from  the  main  program  (Card  B0TR813K) . 
Procrram  Output 

Appendix  B  contains  examples  of  the  full  print  out  and  the 
shorter  print  out  for  the  spherical  wave  reflection  and  a  print  out 
for  a  plane  wave  reflection.  Most  of  the  variables  in  the  output 
are  self-explanatory;  others  which  are  not  so  well  defined  are 
described  below. 

SMALLH  height  h  =*  H  -  d  of  the  charge  above  the  bottom 

DEZFRO  height  d  -  d^  between  the  charge  and  gauge  depths 

D2  reduced  height  dr/R^  from  image  charge  to  gauge 

COSAL  ci  a 

COSTH  cos  9 

SINTH  sin  0 
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DT 

EOT 

T 

STPW 

fi/theta 

PD 

TIME 

PBOT 

PS 

P 

FIMP 


EFLUX 


VMID 

PRE 

RESID 

RFIMP 

REFLUX 

POSIMP 
RPOSIM 
TIME  I 

PTVl 


increment  At  of  the  reduced  time  t 

exp(-At/Gr) 

reduced  time  t 


Vr'‘» 

RiFIlt,/Gr 

incident  pressure  in  psi 

time  in  seconds  relative  to  the  direct  wave  arrival  time 
bottom  reflected  pressure  p^  in  psi 
surface  reflected  pressure  p  in  psi 

total  pressure  p  =  p^  +  Pr  +  pfl  in  psi.  Negative  pressures 

are  cut  off  so  that  P  +  hydrostatic  *  0  • 

total  impulse  I  in  psi-sec  calculated  from  the  equation 
t 


-f 


p  dt,  where  tQ  is  the  minimum  of  6  and  R^/ci 


energy  flux  E  in  in-psi  defined  by  the  equation 


ef  =  (/  |p|  p  at  VG^oeePicj) 


value  of  STPW  at  t  -  At,  R^U  -  At) 
value  of  STPW  ?t  t  -  2 At 


R.  A 

i  / 

reduced  impulse  I/W1  /3 

reduced  energy  flux  Ep/W*  ^3 

impulse  of  the  positive  part  of  the  total  pressure  pulse  p(t) 
reduced  positive  impulse,  POSIMP/W1^3 

time  in  seconds  of  the  PTV,  where  time  is  taken  to  be  zero 
at  the  beginning  of  the  bottom  reflection 
the  PTV  in  ft/sec  induced  by  the  bottom  reflection  in  a 
submerged  target 
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PTV2  the  PTV  in  ft/sec  induced  by  the  bottom  reflection  in  a 
target  at  the  surface 
Time  Independent  FORTRAN  Symbols 


Symbol 

Definition 

Card  Number 

B 

b  =>  Pi/Pa 

B0TR050 

POISR 

Poisson  ratio  o  =  (.5cja-c2)/(c3-c2) 

59 

CSHEAR 

shear  velocity  c*  calculated  from  a 

62 

CSTON 

Stonley  wave  velocity  cgt 

69 

SMALLH 

h  for  zero  slope 

147 

RACTU 

R. 

1 

151 

PH 

negative  of  the  hydrostatic  pressure 

at  depth  dg 

154 

RS 

reduced  surface  reflection  arrival  tine 

158 

W13R 

W1  /Ri 

162 

REDR 

Ri/W1 

163 

THETA 

G 

164 

PACT 

Pp(Ri) 

165 

TACT 

characteristic  time  R^/cj 

166 

THET 

G  in  milliseconds 

167 

A 

bottom  slope  in  radians 

172 

D2ACTU 

dr 

187 

R2ACTU 

Rr 

188 

CTWO 

sin  0cr  =  ci/c8 

213 

R2 

reduced  bottom  reflection  slant  range 

V*! 

214 

THETAR 

2r  -  0(Rr) 

215 

THETR 

Gr  in  milliseco;.is 

216 
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Symbol 

Definition 

Card  Number 

PACTC 

(R  Aj )  G 

r  i# 

n 

B0TR217 

Rl 

(r8aj  G 

s  i' 

218 

THETSR 

G 

8 

219 

SINTH 

sin  0 

221 

COSTH 

cos  8 

222 

D2R2 

cj  6Aa  for  supercritical  reflection 

226 

°»6>*i  for  subcritical  reflection 

228 

SINAL 

sin  6 

cr 

235 

COSAL 

cos  8  ■  cj  a 

cr  1 

238 

SINBE 

sin  8 

crs 

243 

THE 

incident  angle  0  in  degrees 

245 

CR 

plane  wave  reflection  coefficient  K 

261-304 

E2 

phase  shift  <t>  in  radians 

260-307 

EE 

phase  shift  0  in  degrees 

308 

anga 

angle  of  shear  wave  in  bottom  in  degrees 

312, 315 

THONE 

angle  of  compression  wave  in  bottom  in 
degrees 

319,321 

alpha 

0cr  in  degrees 

352 

BETHA 

9 crs  in  de9rees 

358 

SHD2R2 

reduced  arrival  time  of  critically 

refracted  shear  wave 

363-365 

C2 

Cl* 

446 

CBGT2 

cas 

447 

CSHR2 

c4a 

448 

SINTH2 

sin80 

449 

CBSH 

-4c4  3  ,/c8 

450 

42 


NOLTR  71-110 

Symbol 

De finition 

Card  Number 

C2SHR2 

2c4* 

B0TR451 

C4CB 

Cj4b/c2 

452 

Spherical  Wave 

Pressure-Time  Calculations 

463-816 

Symbol 

De finition 

Card  Number 

DT 

increment  At  of  reduced  time 

B0TR476, 500 , 652 

increment  At'  «  At/8 

566,609 

DTI 

original  value  of  At 

478 

DTACT 

2  At/3 

479,501 

EDT 

exp(-At/Gr) 

431,503 

N 

control  parameter  for  pressure  history 

721,804 

VMID 

R.Pr(t  -  At) 

542-662 

STPW 

RiV‘> 

520-690 

PRE 

RiPr(^  "  2i*> 

557-691 

FI 

convolution  integral  Fj 

556-673 

NP 

number  of  subintervals  to  be  used  in 
the  Gaussian  quadrature  integration 
for  Pr(t) 

196,552-693 

V 

v(t)/t  for  t  «  t  -  2a^  where  At7 
c  c 

is  approximately  At/8 

573 

Tl 

t(V) 

578 

T2 

t(.75  V) 

579 

T3 

t(.5  V) 

580 

T4 

t(.25  V) 

581 

U 

u(tc+2  At^)  /t  for  At7  w  At/8 

616 

T2 

t(.25  U) 

623 

T3 

t(.5  U) 

624 

T4 

t ( .75  U) 

62  5 

T5 

t(U) 

626 
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Symbol 

Definition 

Card  Number 

PD 

incident  pressure  p^(t) 

B0TR702 

PS 

surface  reflected  pressure  pa(t) 

704 

PBOT 

bottom  reflected  pressure  pr(t) 

707 

P 

total  pressure  p  ■  p.  +  p  +  p  .  Negative 

i  r  8  ^ 

values  of  p  are  cut  off  at  p  +  hydrostatic  *  0.  713 

Impulse  and  Enerov  Flux  Calculations 

717-767 

De finition 

Card  Number 

XP 

maximum  of  pressure  p  and  zero 

BOTR720 

PMID 

pressure  p  of  even  numbered  time  t-At 

724 

XPMID 

maximum  of  PMID  and  zero 

725 

PPRE 

pressure  p  at  odd  numbered  time  t-2At 

763 

XPPRE 

maximum  of  PPRE  and  zero 

764 

PEND 

pressure  p(t)  at  odd  numbered  time 

757 

XPEND 

maximum  of  PEND  and  zero 

758 

Variables  Used 

in  the  PTV  Calculation  and  in  Plottinq 

Symbol 

Definition 

Card  Number 

XX 

storage  array  for  time  in  microseconds 

for  CALCOMP  plot.  Here  time  is  zero 

at  the  arrival  of  the  direct  wave. 

B0TR800.889 

YY 

storage  array  for  the  total  pressure 

p  for  plot 

801,890 

IPMAX 

number  of  plot  points  stored  in  XX  and 

YY  arrays 

807 

QX 

the  array  in  which  the  time  in  seconds 

is 

stored  for  PTV  calculations.  This  time  is 

802E.813H 

zero  at  the  beginning  of  bottom  reflection,  and  891E 
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Symbol 

QY 


TIMER2 


XT  3 


XT4 


XT  5 


COSA 


Definition  Card  Number 

the  array  in  which  the  bottom  reflection 
pressure  pr(t)  is  stored  for  PTV  calcula¬ 
tions.  If  pr(t)  is  negative,  the  value 
stored  in  QY  is  calculated  so  that 

B0TR802F , 8131 

Pr(t)  +  hydrostatic  *  0  and  891F 

arrival  time  of  the  peak  or  singularity 
of  the  bottom  reflection.  Time  in  this 
case  is  measured  from  the  beginning 
of  the  bottom  reflection.  813C 

signals  the  approach  of  the  bottom 
reflection  singularity.  The  value 
TIMER2  -  2 At  is  used.  The  symbol  T3 
is  used  for  this  variable  in  SUBROUTINE 
PTV.  813D 

The  earliest  time  at  which  the  trans¬ 
lational  velocity  is  to  be  calculated. 

The  symbol  T4  is  used  instead  of  XT4  in 

SUBROUTINE  PTV.  813E 

the  largest  value  of  time  at  which  the 

translational  velocity  is  to  be 

calculated.  The  symbol  T5  is  used  instead 

of  XT 5  in  SUBROUTINE  PTV.  813F 

cosine  of  the  angle  which  the  bottom 

reflection  ray  makes  with  the  water  surface 

or  a  line  parallel  to  the  surface  if  the 

gauge  position  is  below  the  surface  813J 
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Symbol 

Definition 

Card  Number 

PTS 

the  number  of  times  at  which  the  trans¬ 

lational  velocity  is  calculated  in  the 

initial  search  for  the  PTV.  In  the  call 

to  SUBROUTINE  PTV  the  value  PTS  =  30. 

is  used* 

B0TR813L 

Plane  Wave 

Bottom  Reflection  Variables  B0TR819-905 

Symbol 

Definition 

Card  Number 

SW 

direct  wave  response  p^(t)/pp 

B0TR852 

PRFL 

bottom  reflection  response  pr (t) /p£ 

857-877 

TBTH 

(t  -  tc)/Gr 

863 

XE1 

exp ( -TBTH)  E! ( -TBTH) 

865 

XEI 

exp ( -TBTH)  Ei(  TBTH) 

874 

3.2  FORTRAN  Symbols  of  SUBROUTINE  STONL 

Symbol 

De finition 

Y2 

FY 

=  c8st 

Equation  (2.4.2)  which  defines  ya 

CK 

increment  of  ya  =  ya/1000 

CSTON 

Stonley  wave  velocity  c  t 

3.3  FORTRAN  Symbols  of  SUBROUTINE  STPWA 

Symbol 

Definition 

Card  Number 

TR 

c*  Tm 

STPA013 

V 

Ci(cra  -  '>ma)i/s 

14 

Cagniard -Rosenbaum  Method,  CONTR  =  0. 

P(9) 

0.  for  precursor 

20 

XM 

Ci  M 

21 
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Symbol 

Definition 

Card  Number 

P(l) 

Ci  (a  -  M) 

STPA022 

P(2) 

4(cr3  -  T*)1//8sin  9/(a  -  M  ) 

23 

P(5) 

Ci  (a  +  M) 

24 

FACTOR 

2  ^2  b  Rj/tiRj. 

25 

STPW 

Ri 

27,43 

Complex  Arithmetic  Method,  CONTR  =  3,0 

P(l) 

CiX  =  0, 

35 

P(2) 

C1  Tm 

m 

36 

P(3) 

Ci  /c3  =  ci  yi 

37 

P(4) 

ci  y2 

38 

FACTOR 

TTRr/2Ri 

39 

ANS2 

a8 

41 

3.4  FORTRAN  Symbols  of  SUBROUTINE  STFWB 

Symbol 

Definition 

Card  Number 

TR 

C1  Tt« 

TO 

STPB014 

Cagniard-Rosenbaum  Method,  CONTR  =0, 

P(9) 

1.0  for  main  waye 

20 

XK 

ci  K 

m 

23 

XL 

cf  L 

24 

P  (7) 

XK 

26 

P(8) 

XL 

27 

P(ll) 

Ci  2  D 

29 

P  ( 12 ) 

Ci4E 

30 

P(13) 

Ci8F 

31 

FACTOR 

2bRi/TTRr 

32 
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Definition 

Card  Number 

TERMl 

(R.Ar)d  -  b)/(l  +  b) 

STPB038 

SIGM 

cj  a 

42 

STPW 

¥r(t) 

43,47,79,85,100 

XGl 

cigi 

60 

XG3 

Ciga 

61 

XG4 

ci  g4 

62 

XSA 

ciaa/R.8 

63 

XSF 

ci4fAi4 

64 

XNUM 

numerator  of  T 

65 

XDEN 

ci8  times  the  denominator  of  T 

66-67 

RES  ID 

R.  A 

1 

68-69 

TERMl 

Ri(l/Rr  +  a) 

71 

SIG2 

Cl  (<*»  -  cT1)>/* 

78 

Complex  Arithmetic  Method,  CONTR  =»  3.0 

P(l) 

CiX  =  Cx  COS  tt  ( -  cr3)l/2 

93 

P(2) 

Cj  T 

m 

94 

P(3) 

CiYi  =0 

95 

P(4) 

cj  ys  =  Ci  sin  0 

96 

FACTOR 

hr  /2R . 
r  l 

97 

ANS2 

As 

99 

3.5  FORTRAN  Symbols  of  FUNCTION  ONE 

Symbol 

De finition 

Card  Number 

Precursor 

Variables 

X 

integration  variable 

ONE  008 

w 

Cj  u> 

20 
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Symbol 

Definition 

Card  Number 

XC2 

Ci8(a8  -  w2) 

ONE  021,28 

FX 

Fx  defined  by  Equation  (2.2.4) 

23 

Main  Reflection  Variables 

W 

integration  variable  w  =  Ci w 

27 

FX 

Fx  defined  by  Equation  (2.2.21) 

32 

Variables 

for  Rigid  Bottom  Precursors  and  Main  Waves 

FRCS 

c4  8  ( uf*  +  cT2  -  ci*8  ) 

41 

XA 

A 

X 

43 

XB 

Bx 

44 

XC 

Cx 

50 

ONE 

rigid  bottom  integrands  defined  in 

Equations  (2.2.3),  (2.2.14),  and 

(2.2.32)  and  the  integrand  of  the  first 

integral  in  Equation  (2.2.37). 

47,51 

ONE 

fast  non-rigid  bottom  integrands  of 

Equations  (2.2.3)  and  (2.2.21) 

55 

3.6  FORTRAN  Symbols  of  FUNCTION  0NE1 

Symbol 

Da  finition 

Card  Number 

X 

integration  variable  x  =  CiTu 

ONE 1008 

XAB 

Ci E  A  **  Ax 

15 

XBB 

ciBB  *  Bx 

16 

XCB 

ci 6  c  =*  cx 

17 

FAB 

F^  defined  by  Equation  (2.2.38) 

19 

FBB 

F_  defined  by  Equation  (2.2.39) 

D 

20,21 

0NE1 

integrand  F.F^  of  the  second  integral 

A  jO 

in  Equation  (2.2.37) 

49 
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3.7  FORTRAN  Symbols  in  FUNCTION  TWO 

Symbol 

Definition 

Card  Number 

X 

integration  variable  x  =  ci<o 

TWO  007 

FAB 

Fa  defined  by  Equation  (2.2.28) 

13 

FBB 

Fn  defined  by  Equation  (2.2.29) 

D 

14,15 

TWO 

integrand  in  Equation  (2.2.27) 

17 

3.8  FORTRAN  Symbols  in  FUNCTION  SEVEN 

Symbol 

Definition 

Card  Number 

Z 

integration  variable  z  »  Ciy 

SEVN007 

V 

Ci  u 

14 

RCOE 

non-rigid  bottom  K  defined  by  Equation 

(2.3.5)  22 

rigid  bottom  K  defined  by  Equation  (2. 

3.4)  29 

F 

F  defined  in  Equation  (2.3.12) 

34 

SEVEN 

integrand  of  Equation  (2.3.10) 

•9C 

Aa  defined  by  Equation  (2.3.11) 

RT5 

Kj 

39 

U1 

x  +  iyi 

40 

U2 

Cl  3  *1  8 

41 

U3 

Cl  U>1 

42 

XB 

-Ci T  cos  9 
in 

43 

Symbol 

Definition 

Card  Number 

A 

array  a^  in  Equation  (2.5.19) 

EXE1008 

B 

array  b^  in  Equation  (2.5.19) 

9 

C 

array  a^  in  Equation  (2.5.13) 

10,11 

X 

X 

12 
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Symbol 

Definition 

Card  Number 

ANS 

exp(x)  Ex  (x)  for  x  ^  1 

EXE1014, 15 

exp(x)  Ei  (x)  for  0  *  x  <  1 

17,18 

3.10  FORTRAN  Symbols  in  SUBROUTINE  EXEI 

Symbol 

Definition 

Card  Number 

Y 

X 

EXEI006 

A 

array  of  1/nnJ  for  n  =  2,3,.. ,7  in 

Equation  (2.5.20) 

9 

U 

sum  of  the  series  in  Equation  (2.5.20) 

11 

ANS 

exp(-x)Ei(x)  using  Equation  (2.5.20) 

12 

ANSI 

integral  in  Equation  (2.5.21)  evaluated 

using  the  Gaussian  quadrature  of  FUNCTION 

FGI 

15 

ANS 

exp(-x)Ei(x)  using  Equation  (2.5.21) 

16 

3.11  FORTRAN 

Symbols  in  FUNCTION  EXPO 

Symbol 

De finition 

X 

integration  variable  y  in  Equation  (2.5. 

21) 

P(l) 

•x  =  (t  -  Rr/c1)/Gr 

EXPO  integrand  exp(y  -  x)/y  in  Equation  (2.5.21) 

3.12  FORTRAN  Symbols  in  FUNCTION  FGI 
Symbol  Definition 

A  lower  limit  of  integration 

B  upper  limit  of  integration 
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Symbol 

K 


F 

P 

FGI 


Definition 

number  of  subintervals  into  which  the  integration 
interval  (A#B)  is  divided.  The  integral  in  each 
subinterval  is  evaluated  using  a  4  point  Gaussian 
quadrature. 

integrand  of  the  integral  to  be  evaluated 

array  used  to  transfer  parameters  to  the  function  F 

value  of  the  integral  of  F  between  A  and  B 
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APPENDIX  A 

FORTRAN  IV  LISTING  OF  THE  BOTTOM 
REFLECTION  PROGRAM  BOTREF 


B0CH0TR »T600»CM70000. 52400311 t 047. BRITT, 
ATTACH(ABC.nOlBIN) 

COPYN(O.OEF.ABC) 

FTN(L) 

LOAD (L60) 

REQUEST. TAPE99.L0. (CALCOHP/RING) 

OEF. 

'  RECORD  SEPARATOR  =(7-8-9)  PUNCH  IN  COLUMN  1 
REWIND(ABC) 

CALCH1.13.ABC 

'  RECORD  SEPARATOR  =(7-8-9)  PUNCH  IN  COLUMN  1 


C 

C 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


PROGoam  >^0TREF( INPUT. OUTPUT. TAPES* INPUT »TAPE6bOUTPUT.TAPE99)  BOTROol 

B0TR002 

BOTTOM  REFLECTION  PROGRAM  <C0C  6400  COMPUTER)  BOTR003 

OIMEmSION  XX( 1000) *YY (1000)  BOTR004 

COMMON  /OXY/QX(1000) tQY (1000)  B0TR004A 

COMMON  8,COS4L.C0STH,R2,SINBE.SInTH.C*ATER,CB0T.CSHEaR.CsT0N,R£SID  80TR005 

COMMON  C2.CB0T2.CSHR2.CBSH.C2SHR2.C4CB.SINTM2 
AOATF  ■  DATE (O) 

ICASf=1 
PI*3. 1415926 


READ  INPUT  DATA  (FORMAT  —  BF10.S) 

WCH - EXPLOSIVE  CHARGE  WEIGHT  (LBS  OR  Kj) 

CWAtEP— SOUNO  VELOCITY  OF  WATER  (FT/SEC) 

CBOt— --SOUNO  VELOCITY  OF  the  BOTTOM  MATERIAL  (FT /SEC ) 

CSHEAR--IF  CSHEAR  OT  0.5.  IT  IS  THE  SHEAR  WAVE  PROPAGATION 

VELOCITY  OF  THE  BOTTOM  (FT/SEC) ,  IF  CSHEAR  LE  0.5.  IT  IS 
THE  DIMENSIONLESS  POISSON  RATIO. 

RHO  jAt— DENSITY  OF  water  (GM/CC) 

RHOROT— DENSITY  OF  BOTTOM  MATERIAL  (GM/CC) 

PRErOE— COEFFICIENT  OF  PRESSURE  SIMILITUDE  EQUATION  (PSD 
25-.-.. -CONTROL  PARAMETER.  Z5  GREATER  THAN  ZERO  RESULTS  IN  A 
SHORTER  PRINT  OUT. 

PREFXP— EXPONENT  OF  PRESSURE  SIMILITUDE  EQUATION 
THErOE —COEFFICIENT  OF  TIME  CONSTANT  SIMILITUOE  EQUATION  (SEC) 
THEEXP— exponent  OF  time  CONSTANT  SIMILITUDE  equation 
STEOS— NUMBER  OF  POINTS  IN  P-T  CURVE  FOR  ONE  TIME  CONSTANT 
OURAT— DURATION  OF  PRESSURE  TIME  HISTORY  IN  MULTIPLES  OF  THE 
TIME  CONSTANT.  (IF  NEGATIVE.  ITS  ABSOLUTE  VALUE  IS  THE 
OURATION  AFTER  THE  ARRIVAL  OF  THE  BOTTOM  REFLECTION 
PEAK.  IF  POSITIVE.  IT  IS  THE  DURATION  AFTER  THE 
DIRECT  WAVE  ARRIVAL.) 

XI-.— — CALCOMP  PLOT  SCALING  PARAMETER  FOR  THE  Y-AXIS  (PSI  PER 
INCH  OF  GRAPH) 

- - SCALING  PARAMETER  FOR  THE  X-AXIS  (MIcROsECOnOS  PER 

INCH  OF  GRAPH) 

SLOPE— SLOPE  OF  BOTTOM  FROM  CHARGE  TO  GAUGE  (DEGREES) 

AOOlTlCNAL  DATA  IS  READ  IN  STATEMENT  4  (CARO  B0TP0B9) 


BOTR006 

BOTR007 

80TR008 

80TR009 

BOTROlO 

BOTR011 

B0TR012 

B0TR013 

BOTR014 

BOTR015 

80TR016 

BOTR017 

B0TR01B 

0OTRO19 

0OTRO2O 

0OTRO21 

B0TR022 

B0TR023 

B0TR024 

B0TR025 

BOTR026 

B0TR027 

BOTR028 

BOTR029 

80TR030 

B0TR031 

BOTR032 

BOTR033 

B0TR034 

B0TR035 

BOTR036 

B0TR037 

B0TR036 

B0TR039 
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3  READ (5.554) WCH.CwATEr.CSOT »CSHE AR.RHOWAT .RH0B0T .PREC0E.Z5 
READ  (5.554)  PREEXP.THEC0E.THEEXP.STEPS.DURAT.X1 ,X2» SLOPE 
FORMATS  ARE  LISTEO  AT  THE  END  OF  THE  PROGRAM.  CAROS  ROTR907-1042 

DUPAT  IN  RRInT  OUT  IS  THE  DURATION  AFTER  THE  DIRECT  ARRIVAL 
STORF  ORIGINAL  OURAT 

XOURATaOuRAT 

RsRHnWAT/RHOROT 

POISSON  RATIO 

IF  CSHEAR  IS  0.5  FT/SEC  OR  LESS.  THE  POISSON  RATIO  POISR  IS  SET 
EQUAi.  TO  CSHEAR  ANO  THE  SHEAR  VELOCITY  IS  CALCULATED. 

!F(CsHEAR.LE.O.)  GO  TO  39 
IF  (CSHEAR. LE.0.5)  GO  TO  42 

44  POISQ»(0.5»CROT*»2-CSHEAR**2)/(CPOT«*2-CSHEAR«42) 

GO  Tn  41 
42  POISo»CShEAR 

C5HEARaCROT»SORT( (0 .5-PO ISR ) / ( 1. -POISR) ) 

GO  To  41 
39  POlSoaOaH 


STONi EY  "AVE  PROPAGATION  VELOCITY 
41  CALL  STO^L 


THE  GEOMETRY  IS  NOW  READ  IN  (FORMAT  —  8F10.5) 

BIGh— — WATER  DEPTH  AT  THE  CHARGE  (FT).  ALSO  USED  AS  A 
CONTROL  VARIAHLE.  SEE  CARDS  80TR094-97  BELOW 

0— - DEPTH  OF  THE  CHAHGE  BELOW  THE  WATER  SURFACE  (FT) 

OGAil— — OEPTH  OF  THE  GAUGE  (FT) 

SMAlLR— HORIZONTAL  RANGE  BETWEEN  CHARGE  AnD  GAUGE  (FT) 

THOvAL— DESIRED  RATIO  BETWEEN  THE  BOTTOM  REFLECTION  INCIDENT  AND 
CRITICAL  ANGLES.  (IF  THOVAL  LE  0.*  THE  INPUT  GEOMETRY  IS 
NOT  CHANGED.)  SEE  APPENDIX  C  OF  NOLTR  71-110. 

Zl— — — PARAMETER  which  SELECTS  THEORY.  (SEE  CARDS  B0TR391-394) 

Z2 - ARRIVAL  TIME  OF  THE  MAIN  BOTTOM  REFLECTION  AFTER  THE 

DIRECT  ARRIVAL  (MICROSECONDS) .  GEOMETRY  IS  UNCHANGED  IF 
Z2  LE  0,  SEE  APPENDIX  C  OF  NOLTR  71-110. 

23— ——PLOT  CONTROL  PARAMETER,  A  CALCOMP  PLOT  TAPE  IS  GENERATED 
IF  Z3  IS  ZERO  . 

4  READ (5,554) 8 IGH.D.DGAU. SMALLR, THOVAL. Z1.Z2.Z3 

FORMATS  ARE  LISTEO  AT  the  ENO  OF  THE  PROGRAM,  CAROS  ROTR9O7-1042 

aftep  completion  of  each  case  program  control  returns  to 
STATFMENT  4,  oepenoing  on  bigh,  the  CALCULATIONS  are  continued  as 
FOLLOWS 

IF  BIGH  «  0  PROGRAM  STOPS 

IF  BIGH  IS  POSITIVE  COMPUTATION  CONTINUES  USING  THE  PRESENT  INPUT 
IF  BIGH  Is  NEGATIVE  PROGRAM  TRANSFERS  TO  STATEMENT  3  WHERE 
ANOTHER  SET  OF  CHARGE,  PHYSICAL  CONSTANTS.  ETC.  ARE  READ, 

IF (BIGH) 3,1000.6 
1000  STOP 

6  w«ITc(6.510) AOATE 

J  FORMATS  ARE  LISTEO  AT  THE  ENO  OF  THE  PROGRAM,  CAROS  BOTR907-1042 
:  ADDITIONAL  OATA  IS  READ  IN  FOR  PTV  CALCULATION  (FORMAT  —  8F10.5) 


BOTR040 

BOTR041 

80TR042 

BOTR043 

BOTR044 

80TR045 

BOTR046 

80TR047 

BOTR04B 

BOTR049 

BOTR050 

B0TR051 

B0TR052 

80TR053 

BOTR054 

80TR05S 

B0TR056 

80TR057 

BOTRO58 

BOTR059 

80TR060 

BOTR061 

80TR062 

BOTR063 

80TR064 

80TR065 

80TR06P 

BOTROG7 

e0TR068 

BOTR069 

BOTR070 

80TR071 

B0TR&72 

80TR073 

BOTR074 

BOTR075 

B0TR076 

BOTR077 

B0TR078 

BOTR079 

BOTRO0O 

BOTROfll 

80TR082 

BOTROB3 

BOTROB4 

BOTRO05 

BOTRO06 

BOTRO07 

BOTRO08 

BOTRO09 

BOTR090 

BOTRO9I 

BOTR092 

B0TR093 

B0TRC94 

BOTR095 

BOTR096 

BDTR097 

BOTR09B 

BOTR099 

BOTR100 

BOTRlol 

BOTRIOIA 

BOTRIOIB 
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C 


RAOlUS— CYLINDER  RA0IUS  IN  FEET.  THIS  IS  THE  DRAFT  OR  C&0SS- 
SECTIONAL  RADIOS  OF  THE  TARGET  VESSEL. 

APRtNT— controls  PRINTING  in  subroutine  ptv.  the  ITERATIONS  to 
OBTAIN  THE  PTV  ARE  PRINTEO  OUT  if  APRlNT  .LE.  0. 

REAOfS.S'U)  RADIUS* APRINT 

IPTV.O 

AaO. 

WRITF (6,550) ICASE 
WRITF(6,511)BIGH 
WRITf (6,513) 0 
WRITF(6,M2'0GAU 

FORMATS  ARE  !ISTE0  AT  the  END  OF  THE  PROGRAM,  CAROS  ROTR907-1042 
CHANGE  OF  DEPTH  OF  EXPLOSION  ANO  OF  GAUGE  FOR  GIVEN  ARRIVAL 
TIME  22  OF  THE  BOTTOM  REFLECTION  PEAK  AFTER  THE  DIRECT 
ARRIVAL.  THE  BOTTOM  REFLECTION  INCIDENT  ANGLE  ANO  SmALLR 
ARE  UNCHANGEO.  SEE  APPENOIX  C  OF  NOLTR  71-110, 

GEOMETRY  IS  UNCHANGEO  IF  Z2  IS  NEGATIVE  OR  ZERO. 


-D2ACTU-  IS  THE  ORIGINAL  TOTAL  DISTANCE  BETWEEN  THE  GAGE  AND  THE 
IMAGf  CHARGE,  -R2ACTU-  IS  THE  ORIGINAL  TOTAL  SLANT  DISTANCE  FROM 
THE  CAGE  TO  THE  IMAGE  CHARGE  AND  -HG-  IS  THE  NEW  VALUE  OF  THE 
HEIGHT  OF  THE  GAGE  ABOVE  THE  BOTTOM.  NEW  VALUES  FOR  -0-  AND 
-OGAil-  ARE  CALCULATED. 

IF ( Zp.LE.O. >  GO  TO  2 

1  02ACtU*2.*<BIGH-0).0-0GAU 
R2ACxU«Sf}RT(SMALLR**2.02ACTU»«2) 

0ELR-Z2*  ( 1  .E-06)  *'CWATER 
C12002«0ELR*(2.«R2ACTU-0ELR)/02ACTU»*2 
HG«0 .5«02ACTU* ( 1 .-SORT  1 1 .-0R2002 ) ) 

OGAUaBIGH-HG 

0»8IgH-02ACTU*HG 
WRITr (6,553)  Z2 
WRITF (6,513)0 
WRITf(6,512)0GAU 

change  of  geometry  to  obtain  the  desired  ratio 

BETWEEN  INCIDENT  ANO  CRITICAL  ANGLE-THOVAL 

GEOMfTRY  IS  UNCHANGEO  IF  THOVAL  IS  LESS  THAN  OR  EQUAL  TO  ZERO. 
SEE  APPENOIX  C  OF  NOLTR  71*110. 

2  IF(ThOVAL.LE.O.)  GO  TO  5 

7  THaTHOVAL*  ASIN  ( CWATER/CBOT) 

02ACTU»SmALLR*C0S ( TH) /SIN (TH) 

0a2.*BIGH-DGAU-02ACTU 
IF(BTGH-n)  8*9,9 
B  D  a  PIGH 

DGAU  a  HIGH  -  02ACTU 
9  WRITf (6,537) 

WRITE (6 ,512) DGAU 
wRITf (6,513)0 
GEOMfTRY 
5  SMALLHaBIGH-D 

-RACTU-  IS  the  SLANT  DISTANCE  BETWEEN  charge  AND  GAG?;. 

RACTu  a  SQRT< (0*D3AU)«#2*SMALLR**2) 

calculate  hydrostatic  pressure  -PH» 

PHa-i4.7«DGAlJ/33.0-14.7 


BOTRlOlC 

botrioid 

BOTRIOiE 
BOTRlOlF 
BOTRlOiG 
BOTR101H 
BOTRlOlI 
BOTR101J 
BOTR101K 
BOTR1 o2 
8OTR103 

botrio* 

BOTR1o5 

BOTR106 

BOTR107 

BOTR108 

BOTR109 

botruo 

BOTRlil 
BOTRU2 
BOTR1 l3 
80TRU4 
B0TR115 
80TRU6 
B0TRU7 
BOTRllB 
BOTR 1 i 9 
80TR120 
BOTR121 
80TR122 
80TR123 
B0TR124 
BOTR125 
B0TR126 
BOTH  127 
90TR1PB 
B0TR129 
BOTR130 
BOTR131 
B0TR132 
BOTR133 
BOTR134 
BOTR135 
80TR136 
BOTR137 
BOTR130 
B0TR139 
BOTRUO 

botrki 
B0TR142 
BOTRH3 
BOTRH4 
BOTR 143 
B0IR146 
B0TR147 
60TRU8 
BCTRU9 
BOTR150 
B0TR151 
80TR1P2 
BOTR 153 
B0TR15* 
B0TR155 


A-3 
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-RS-  is  the  reduced  arrival  time  of  acoustic  surface  reflection, 

HS«S0RT (SMALLR##2* (D^DGAU) #*2 > /RACTU 

EXPONENTIAL  PULSE  PEAK  PRESSURE  AND  TIME  CONSTANT  CALCULATED 

W13R«WCH«»(l./3.) /RACTU 
REDR«1./W13R 

THETA=THEC0E«(V(13R>««(1,4THEEXP>»CWATER 
PACT*PREC0E*(W13R)»»PREEXP 
TACT«RACTU/CWATER 
THET«THETA»TACT*1000. 

IF  (SLOPE  .EG.  0.)  GO  TO  10005 

CHANGE  OF  GEOMETRY  FOR  SLOPING  BOTTOM 

A  «  SLOPE/57.29578 
HG  «  BIGH-DGAU 
HI  •  SMALLH«COS(A) 

H2  »  (HG-SMALLR«TAN(A) >»COS(A) 

IF (  (H2.LT.D.0 ) cOR. (SMALLR#TAN (A).GT.BIGH) )  WRITE (6.555) 

WRITE  (6.514)  SMALLR 
WRITE  (6.574) 

SMALLR«SMALLR*COS(A)*(D-DGAU)*SIN(A) 

SMALLH  s  HI 

IF  (H2  ,GT.  HI)  D*D*H2-H1 
DGAU  *  D.H1-H2 
BIGH  «  D*H1 
WRITE  (6.511)  BIGH 
WRITE  (6.512)  DGAU 
DEZERO  =  D-DGAU 
D2ACTUs2.*SMAL!-H*DEZER0 
R2ACTU*SQRT (D2ACTU##2*SHALLR##2) 

INITIALIZATIONS 

1  040  FI»0, 

VMID.O. 

PPE».T, 

IP*I 
NP«4 
Z2DTa4. 

RES I 1*0.0 
EFLU*«0. 

FIKP*0. 

POSI"PbO. 

I?RE<;»^1 

FRREs6. 

XPPRr»0. 

PMIOcO. 

XPMIn»C, 

PO»Q. 

PSsO, 

P3CT.0, 

e*SIO  CONSTANTS  OF  GROUND  WAVE 

CTWOaCw  ATER/C0OT 
R2»R?ACT'J/RACTU 
THET  «R»THETA/R2«*THEEXP 
THETo*THET/R2»*THEEXP 
PACTr«R2*«PREExP/R2 
fil«R<;#«PHEEX» 

ThE7<R*ThETA/R1«*TH£EXP 


B0TR156 
B0TR157 
B0TR158 
BOTR15R 
BOTR160 
B0TR161 
BOTR162 
BOTR163 
B0TR164 
B0TR165 
B0TR166 
B0TR167 
B0TR168 
BOTR169 
BOTR170 
B0TR171 
BOTR172 
BOTR173 
BOTR174 
B0TR175 
B0TR176 
B0TR1T7 
BOTR1 70 
B0TR179 
BOTR180 
BOTR181 
BOTR182 
BOTR183 
B0TR184 
B0TR185 
B0TR186 
B0TR187 
BOTR188 
BOTR189 
BOTR190 
B0TR191 
B0TR192 
90TR193 
B0TR194 
BOTR195 
B0TR1P6 
BOTR197 
B0TR198 
B0TR199 
B0TR200 
BOTR201 
BOTR202 
BOTR203 
BOTR204 
EOTR205 
BOTR208 
BOTR207 
B0TR208 
BOTR209 
B0TR210 
B0TR2U 
BOTR212 
B0TR213 
B0TR214 
B0TR215 
BOTR216 
B0TR21 7 
B0TR218 
B0TR219 
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02»0?ACTH/HACTU 
S1NTh»SMalLR/ (RACTU*R2) 

COSTh»02/R2 
C0SRm»C0STH/R2 
IF(CtWO.GE.SINTH)  GO  TO  2000 
GAM»«50BTn.-CTx0«*2) 

D2R2«R2»fCTWO*SlNTH*GAM»COSTH) 

GO  T n  2005 
2000  02R2»R2 
C 

C  CALCULATION  OF  NEW  OURATlON  IF  OtjRAT  IS  READ  IN  NEGATIVE, 

C  (DUR/.T  IS  EXPLAINED  IN  CAROS  80TR028-32.) 

2005  IF(XngRAT,LT.O.)  OURAT«(R2-l.)/TH£TA-XOURAT 
TSTOo»l,*OURAT*THETA 

C 

SINAl«CWATER/CROT 

SIN2aL»SINAL#*2 

IF1SIN2AL-1. >811* 811*812 

811  C0SAL»SQRT(1»-SIN2AL) 

GO  TO  813 

812  COSALs-O, 

813  SIN2TH»SINTH«M>2 

IF (CSHEAR) 15 » 15  *  14 

14  SlNBE*CWATER/CSHEAR 
SIN2BE»SlNBE**2 

15  THE«57,2958»  ASIN(SINTH) 

C 

C  CALCULATION  OF  PLANE  WAVE  REFLECTION  COEFFICIENT.  PHASE  SHIFT. 

C  AND  ANGLE  OF  S-WAVE. 

C 

IF (CSHEAR)30. 30.50 

C  CALCULATION  FOR  BOTTOM  WITH  NO  SHEAR  STRENGTH  (NON-RIGID) 

C 

30  IF(S.N2TH-SIN2AL)33.32.31 

c  SUPERCRITICAL  REFLECTION  (ANGLE  OF  INCIDENCE  GREATER  THAN  THE 
C  CRITICAL  ANGLE) 

C 

31  E»ATAN(B*SQRT(SIN2TH-SIN2AL)/C0STH) 

CR»1, 

IICAsl 
GO  TO  88 

32  E2«0. 

CR«1, 

I ICA»2 
GO  TO  89 

C  SUBCRITICAL  REFLECTION  (ANGLE  OF  INCIDENCE  LESS  THAN  THE 

C  CRITICAL  ANGLE) 

C 

33  AT»COSTH/ (SQRT (SIN2AL"SIN2TH) #B ) 

CR»(AT-1.)/(AT*1.) 

IICA«3 
GO  TO  89 

C  CALCULATION  FOR  BOTTOM  WITH  SHEAR  (RIGID) 

C  50  CA«C0STH*(1.-2#»SIN2TH/SIN2BE)*»2/B 

CB*4.»C0STH*SIN2TH*(SIN2BE-SIN2TH)/B/SIN2BE*»2 
IF (SIN2TH-SIN2AL) 60.32.51 
51  ATA»CA/SQRT(SIN2TH*SIN2AL> 


B0TR220 

B0TR221 

B0TR222 

B0TR223 

B0TR224 

B0TR225 

B0TR226 

B0TR227 

B0TR228 

B0TR229 

B0TR230 

BOTR231 

B0TR232 

B0TR233 

B0TR234 

BOTR235 

B0TR236 

BOTR237 

BOTR238 

B0TR239 

BOTR240 

B0TR241 

BOTR242 

BOTR243 

B0TR244 

B0TR245 

B0TR246 

B0TR247 

B0TR248 

B0TR249 

BOTR250 

BOTR251 

BOTR252 

B0TR253 

BOTR254 

B0TR255 

B0TR256 

BOTR257 

BOTR258 

B0TR259 

BOTR260 

BOTR261 

BOTR262 

BOTR263 

B0TR264 

BOTR265 

B0TR266 

BOTR267 

BOTP268 

BOTR269 

BOTR270 

BOTR271 

BOTR272 

BOTR273 

80TR274 

BOTR275 

BOTR276 

BOTR277 

BOTR278 

BOTR279 

B0TR280 

90TR281 

B0TR282 

B0TR283 
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IF(SINTH-SIN0E)52*55*57  B0TR284 

52  0T«CB/SQRT(SIn2B£-SIN2TH)  0OTR285 

CR*S0RT { ( ATA»*2* (BT-1 , ) **2) / (ATA*»2* (BT*l , ) *#2) >  0OTR206 


EA«ATAN((1.-0T)/ATA) 

E0"ATAN ( ( 1 «  *0T ) /ATA) 

E2»EA*E0 
IICA«4 
GO  TO  89 
55  0TA«O. 

GO  TO  58 
C 

57  BTA«C0/SQRT(SIN2TH-SIN2SE) 

58  E«ATAN(1,/(ATA*BTA) ) 

ca-1, 

IICA-6 
GO  TO  88 
C 

60  E«0» 

AT-CA/SORT (SIN2AL-S1N2TH) 

0T-C0/SQRT (SIN28E-SIN2TH1 
CR«(AT*BT.l.)/(AT40T*l.) 

1  ICA»7 
C 

88  E2«2 ,*E 

89  EE*57.3959#E2 

IF  (CsHE AR.LE.O, >  GO  TO  92 

90  IF(STNTH-SIN9E)91. 91*92 

91  GAMMas  ASINISINTH/SIN8E) 

ANGAS5? ,2956*0 AMM A 

GO  Tn  95 
C 

92  ANGA.-O. 

C 

c  anglf  of  p-wave 

95  IF(S1NTH-SINAL1293. 293,294 

293  TH0Nf«57.2958*  ASIN (SINTH7SINAL) 

GO  To  295 

294  THONF—O, 

C 

C  FORMATS  are  LISTED  AT  THE  END  OF  THE 

295  WRITF (6*514) SMALLR 
WRITf (6*504 ) WCH 
WRITf (6  *505) CWATER 
WRITf (6.506) CBOT 
WRITf  ',6*546)CSHEAR 
WRITf (6«507)RHOWAT 
WRITF (6 *508) RHO0OT 
WRITf(6,515)PREC0E 
WRIT? (6*503) 75 
wRITr (6,516) PREEXP 
wRITf(6,s17)THECOE 
WRITf (6*518) THEEXP 
wRITr (6,519) STEPS 
WRITf(6,509) OUR AT 
WRITf (6*538) THOVAL 
WRITF (6*504)  XI 
WRITf  (6,585)  X2 
WRITf  (6*500)  SLOPE 
WRITf (6*586)  Z1 
WRITF (6  *587)  Z2 
WRITf (6  *580 ) 73 
WRITf (6*568)  RADIOS 
wRITr (6,549)  APRINT 
WRITF (6*520) 


0OTR287 
BOTR208 
BOTR209 
B0TR290 
0OTR291 
B0TR292 
B0TR293 
B0TR294 
B0TR295 
0OTR296 
801R297 
B0TR298 
B0TR299 
0OTR3OO 
0OTR3O1 
0OTR3O2 
BOTR303 
BOTR304 
B0TR305 
0OTR3O6 
BOTR307 
B0TR308 
B0TR3Q9 
BOTR310 
0OTR3H 
B0TR312 
BOTR313 
B0TR314 
B0TR315 
B0TR316 
0OTR317 
B0TR318 
80TR319 
BOTR320 
0OTR321 
80TR322 

PROGRAM,  CARDS  BOTR907-1042  BOTR323 

0OTR324 

B0TR325 

B0TR326 

BOTR327 

B0TR328 

80TR329 

BOTR330 

0OTR331 

0OTR332 

0OTR333 

0OTR334 

0OTR335 

BOTR336 

0OTR337 

80TR33B 

0OTR339 

BOTR340 

B0TR341 

B0TR342 

P0TR343 

B0TR34* 

B0TR344A 

B0TR344B 

B0TR345 
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WRITE (6*521 ) THE 
WRITE (6»573)CST0N 
WRITE (6 *545 5 POl SR 
WRITE  (6.547)  RS 

IF(I,-SINAU17.16.16 

16  ALPHA*57»2958*  ASIN(SINAL) 
write (6*522) alpha 
GO  TO  18 

17  WRITE (6*541 ) 

18  IF  (CSHEAR)49. 49*45 

45  IF  ( 1 ,-SINBE) 47*46*46 

46  0ETHAs57.2958«  ASIN(SINBE) 

WRITE (6.542) 8ETH A 

ARRIVAL  TIME  OF  CRITICALLY  REFRACTEO  SHEAR  WAVE 
IF(STNTH,LT.SIN8E)  -SH02R2»0. 

IF(STNTH.6E*SINBE  )  SH02H2«(SMALLR*SINBE*02ACTU#SQRT(1.-SIN2BE) ) 
1  /RArn 

wRITf (6 .579)  SHD2R2 

GO  Tn  49 

47  WRI Tf ( 6*543 ) 

49  WRITF (6.597JTH0NE 

wRl Tf (6  *592) CR 
wRITf(6»594) ANGA 
WRITf  (6.S9DEE 
WRITf (6»k23)02R2 
wRITf (6  *533) R2 
wRITf (6 *525) R AC TU 
wRITf (6 .^02 )  REOR 
WRITF(6,R26)TACT 
WRJTf(6*5?7)PACT 
WRITf (6,528) THETA 
WRITf (6*539) T HE T 
WRITf (6*548)  THETAR 
wRITf (4  * 549 )  ThETH 
WRITf (6,535) 

WRITf (6,551  )SMALLH*0E7FR0*02*C0SAL*C0STH.SINTh 
wRITf (6,532) 

FORMATS  are  LISTED  AT  THE  END  OF  THE  PROGRAM,  CArDS  BOTR907-1042 

SELECT  tON  OF  THE  THEORY 

Z1*0.  ROSENBAUM  METHOD 
Zl*l.  PLANE  WAVE  APPROXIMATION 
21*2.  NOT  USED  In  THE  PRESENT  PROGRAM 
71*3.  COMPLEX  ARITHMETIC  METHOD 

Z l*AaS ( Z 1 ) 

IF(Zl-l,)R00, 801,802 
802  IF  (Zl-3, ) 803.804,805 

SPHERICAL  WAVE  CAGMARO-RCSENBAUM  THEORIES 

800  WRITf (6 ,567 ) 

IE (CsH£A9) 82 0, R20.821 

C 

C  NON-uIGIi'  bottom 

820  IF(STNAL-l. >830*831*832 
C 

C  FAST  8CTTOM 

830  WRITf (6*560) 


B0TR346 
B0TR3*7 
B0TR348 
B0TR349 
BOTR350 
B0TR351 
BOTR352 
B0TR353 
B0TR354 
B0TR355 
B0TR356 
B0TR357 
B0TR358 
B0TR359 
BOTR3&0 
B0TR361 
BOTR362 
B0TR3*>3 
B0TR36A 
BOTR365 
B0TR366 
B0TR367 
B0TR368 
B0TR369 
BOTR370 
B0TR371 
B0TR372 
B0TR373 
B0TR37* 
B0TR375 
BOTR376 
BOTR377 
B0TR378 
30TR379 
B0TR3B0 
B0TR3B1 
BOTR382 
B0TR383 
B0TR3B4 
BO'i  R3fl5 
B0TR3B6 
B0TR3B7 
B0TR3B8 
B0TR3B9 
BOTR390 
B0TR391 
B0TR392 
B0TR3  )3 
B0TR39* 
B0TR395 
B0TR396 
B0TR397 

B0TR398 
B0TB399 
BOTR4nO 
B0TR40 1 
B0TR402 
0OTR4O3 
B0TR404 
B0TR4nb 
B0TR406 
80TR407 
BOTR408 
0OTR4n9 
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GO  TO  11 
C 

C  SLOW  BOTTOM 

®32  WRITE (6*561) 

GO  TO  11 
C 

C  NO  REFLECTION 

?3l  WHITE (6*599) 

GO  TO  4 

C 

C  RIGIO  BOTTOM 

82l  IF<SIN8E-1. >841*841. 840 

C 

C  FAST  SHEARWAVE 
§*1  WRITE (6*562) 

GO  TO  11 
C 

C  SLOW  SHEARWAVE 
840  WRITE (6*563) 

GO  TO  11 
C 

c  plane  WAVE  APPROXIMATION 
®01  WRITE (6*565) 

GO  TO  998 
C 

c  Z1  a  2.  IS  NOT  NEEDED  FOR  THE  PRESENT  PROGRAM 

*03  wR I Tr ( 6 * S99 ) 

GO  To  4 
C 
C 

c  COMPI  EX  ARITHMETIC  METHOD 

c 

804  wR;Tm6,S66) 

c 

c  CONSTANTS  FOR  SUBROUTINE  SEVEN 

C2«CvATE»»*2 
CB0T-5«CROT**2 
CSHR?aC$MEAR*«2 
SINTM2«SfNTH*#2 

CBSHc-4.*CSHEAP#»3/C80T 

C2SHn2»2,»cSHEAR*#2 

C4CBsC?**2/C80T*B 

GO  To  11 
C 

C  ?1”*.  IS  NOT  NEEOEO  FOR  PRESENT  PROGRAM 

805  wRITr (6.S99) 

GO  Tn  4 


SPHERICAL  WAVE  CAGM ARO-RCSENBAUM  PRESSURE-TIME  CALCULATIONS 


C 

c 

c 

c 


PHASrS  A'jQ  TIME  STEPS 

11  IF  (0?R2-£>2)  10*20*20 

ANGLP  CF  INCIOENCE  GREATER  THAN  CRITICAL 
10  H«(R’-C2R2)*STEPS/THETA/4, 

IF  <  <r'2B2-l  .0)  .GT.O.  >  Ma  (R2-1.0)»STEPS/THETA/4, 


BOTR410 

B0TR411 

B0TR412 

B0TR4 13 

B0TR414 

BOTR415 

B0TR416 

B0TR417 

B0TR418 

BOTR419 

BOTR420 

B0TR421 

B0TR422 

B0TR423 

BOTR424 

BOTR425 

BOTR426 

B0TR427 

B0TR428 

B0TR429 

BOTR430 

B0TR431 

B0TR432 

BOTR433 

B0TR434 

BOTR43S 

BOTR436 

BOTR437 

B0TR438 

BOTR439 

BOTR440 

B0TR441 

B0TR442 

B0TR443 

B0TR444 

60TR445 

B0TR446 

B0TR447 

BOTR448 

B0TR449 

BOTR450 

B0TR451 

B0TR452 

B0TR453 

B0TR454 

B0TR455 

B0TR456 

B0TR457 

B0TR458 

BOTR459 

BOTR460 

B0TR4f,l 

E0TR462 

B0TR463 

B0TR464 

B0TR465 

B0TR466 

B0TR467 

B0TR468 

B0TR469 

B0TR470 

BOTR471 

B0TR472 

B0TR473 
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M»4*M*5  ,.  _  _  _ 

CALCULATE  DT,  INCREMENT  OF  REDUCED  TIME  T 
12  DT ■  (R2-D2R2 ) /FLOAT (M) /2.  „A,,U. 

IF  (  (D2R2-1.0)  ,GT.0»  )  DT«  <R?-1  .0}  /FlOAT  (M)  /2. 

DT1»DT 

DTACT«2.#DT*TACT/3. 

DST«DT/3. 

EOT«EXP(-DT/THETAR) 

WRITE (6,536) DT  » EOT 

formats' arI’ LXSTEO  AT  THE  END  OF  THE  PROGRAM,  CARDS  JOTWOT-l 042 
7c„  T«  THE  PRINTOUT  CONTROL  PARAMETER.  IF  -25-  GREATER  THAN 
2ERO  ['SHORTER  PRINTOUT  RESULTS,  IF  -25-  EQUALS  ?ERO  THE  NORMAL. 
LONGrR  PRINTOUT  IS  GENERATED, 

IF  (7S.GT,  0.0)  GO  TO  103 
WRITE  (6.530) 

GO  T*  100 
103  wRITc(ft,'-l)l) 

GO  T->  lO'i 

ANGLc  CF  INCIDENCE  LESS  THAN  OR  FQUAL  TO  CRITICAL 

?0  MM»  ( r.2- 1  •  )  *STEPS/THET  A 
mm=2*mr*a 

calculate  dt.  increment  of  reouceo  time  T 
ot» (^2-1 . ) /Float  <mm) 

OTACr«2,*OT»TACT/3. 

0ST»nT/3. 

EDT«FXP(-DT/THETAR) 

WRIT? (6,536) DT, EOT 
wMITf (6 ,532 ) 

IF  (25. 6T, 0.0)  GO  TO  104 
wRITf (6 ,530) 

GO  To  7 OC 
104  WRITf ( 6 ,50 1 ) 

GO  T"  700 


ANGLP  CF  INCIDENT  WAVE  LARGER  THAN  CRITICAL 
R2  LARGER  THAN  D2«2 

ICO  IF  (D?R2-r’,9999J  101. 102, 102 

PRECURSOR  ARRIVES  BEFORE  DIRECT  WAVE 

101  T«D2R2 
STPWsO. 

Ns  1 0 

GO  TO  72 


C 

C 

C 


PRESSURE  CALCULATION  IF  DIRECT  WAVE  ARRIVES  BEFORE  PRECURSOR 


102  T«1 .0 
STPW«0. 

Ns  1 

GO  TO  71 
110  Nsl2 

114  T»T*2.*0T 

IF  (T.LT.D2R2)  GO  TO  71 
117  TsD2R2 

WRITE (6.53A) 

Ns  1 1 

GO  TO  71 


B0TRA74 

B0TR475 

BOTR476 

B0TR477 

B0TR478 

B0TR479 

BOTRA60 

B0TR481 

B0TR4B2 

80TR483 

B&TR4B* 

B0TR485 

B0TR4R6 

B0TR487 

30TR488 

B0TR4B9 

BOTR490 

B0TR491 

B0TR492 

B0TR493 

B0TR494 

B0TR495 

B0TR496 

B0TR497 

B0TR498 

B0TR499 

BOTR500 

BOTRbOl 

BOTR502 

BOTR503 

BOTP50* 

BOTR505 

BOTR506 

BOTR507 

BOTR5O0 

B0TP509 

BOTR510 

B0TR511 

BOTR512 

B0TR513 

B0TR51* 

B0TR515 

B0TR516 

BOTR517 

BOTR51B 

BOTR519 

BOTR520 

80TR521 

BOTR522 

BOTR523 

BOTR524 

BOTR525 

BOTR526 

BOTR527 

BOTR52S 

BOTR529 

BOTR530 

BOTR531 

BOTR532 

BOTR533 

BOTR534 

BOTR535 

BOTR536 

BOTR537 
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C  CALCULATION  OF  THE  PRECURSOR 
C 

150  N«2 
152  T»T*0T 

CALL  STPWA(T.VMlO.Zl.NP) 

C  EVERY  FIFTH  STEP  RECALCULATE  STPW  WITH  TWICE  THE  INTEGRATION 
C  POINTS. NP.  TO  CHECK  THE  INTEGRATION  ERROR. 

IF (MOO ( IP *5) <NE.O )  GO  TO  155 
IF (NP.EQ. 16)  GO  TO  155 
NP2»7«NP 

CALL  STPWA (T  » VMI02*  Z1 *NP2) 

ERR0R«ASS{ (VMID-VMI02)/VMIo2) 

IFJERROR.LT. 0.005)  GO  To  155 
C  DOUBLE  NUMBER  OF  INTEGRATION  POINTS 
NP«NP2 
VMID«VHID2 
155  T*T.0T 

CALL  STPWA(T.STPW.Zl.NP) 

FI«FI«EDT<n*2M  (PRE*EDT*4.«vMID)*EDT*STPW)  *DST 
PREaSTPW 
GO  TO  70 

159  IF ( (IPRES.LT.O) .OR. (T.LT. <R2-6.1*DT> ) )  GO  TO  150 
CALCULATION  of  precursor  near  SINGULARITY 

200  DUOT/ZZDT 
M«(R2-T)/DT/4. 

M*4#H #5 

DT« (R2-T) /FLOAT (M) /2, 

OTACT«2.*DT# TaCT/3. 
dst-ot/3. 

EOT«EXP(-DT/THETAR) 

N«9 

201  IF( (IPRES.LT.O) .OR. (T.LT. (92-3. 1*DT) ) )  GO  TO  152 
TR1.T/R2 

V-SQRT(1.-TR1««2) 

0SV»V/12.*R2 
TR2»SQRT  ( 1 ,  -  (0 .75*V )  *«*2) 

TR3»SQRT(1.-(0.5*V)*«2) 

TR4»SQRT ( 1 ,- ( 0 .25*V) ••Z) 

Tl-T 

T2«R2*TR2 
T3«R2*TR3 
T*«R2*TR4 

EDTUEBP(-(T3-T1)/THETAR) 

EDT2*EXP  <• (T3-T2) /THETAp) 

EDT3«EXP(-(R2-T3)/THETAR) 

EDTA«EXP(-(R2-TA)/THETAR) 

202  CALL  STPWA (T2.VMI0.Z1. 16) 

CALL  STPWA (T3.STPW.Z1. 16) 

FUFUEDT1  ♦  (PRE*E0T1*V/TR1  ♦3,«VMIo*EDT2*V/TR2.STPW*0 .5*V/TR3 )  *DSV 
C 

PRE«S7PW 

T«T3 

N»3 

GO  TO  70 
C 

210  CALL  STPWA (T4.VMID.Z1 » 16) 

FUFI*E0T3* (PRE*EDT3*0 .5*V/TR3*VMI0#EdTA*V/TR4) «0SV 
PRE"0. 

T»R2 

STPW»CR* (1.E*30)*SIGN(1,»£E> 

WRITE (6.540) 


B0TR538 

B0TR539 

BOTR540 

80TR541 

B0TR542 

B0TR543 

R0TR544 

B0TRS45 

B0TRB46 

B0TR547 

B0TR548 

80TR549 

B0TR550 

B0TR551 

BOTR552 

B0TR553 

B0TR554 

B0TR555 

B0TR556 

B0TR55? 

B0TR558 

B0TR559 

BOTR560 

B0TR561 

B0TR562 

BOTR563 

R0TR564 

B0TR565 

B0TR566 

50TR567 

B0TR568 

B0TR569 

GOTR5TO 

BOTR571 

B0TR572 

B0TR573 

B0TR574 

B0TR575 

B0TR576 

B0TR577 

B0TR578 

B0TR579 

BOTR580 

B0TR581 

BO. R582 

B0T&.I81 

B0TR5B4 

B0TR585 

B0TR586 

B0T8587 

B0TR588 

B0TR589 

BOTR590 

80TR591 

B0TR592 

B0TR593 

B0TR594 

BOTR595 

80TR596 

BOTR597 

B0TR598 

B0TR599 

BOTR600 

BOTR601 
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N»4 

60  TO  70 


CALCULATION  OF  MAIN  BOTTOM  REFLECTION  NEAR  SINGULARITY 

300  IF(DT.GT.(0Tl/ZZ0T/2.) )  GO  TO  301 
DT*DTl/ZZDT/2 • 

0TACT=2.*DT*TACT/3. 

DSTsDT/3. 

EDT-EXP (-0T/THETAR) 

NP*16 

301  T6=R2*4.*DT1 
DTR*DT/R2 

U-SQRT ( ( 1 .♦2,*DTR>  **2-1 . ) 

DSU«U/12.*R2 

TR2=SORT(1.*(0.25*U)**2) 

TR3xSORT(1.»(O.SO«U}**2) 

TR4sS0RT ( 1 . ♦ ( 0 .7S*U) »#2 ) 

TR5=1.+2.*DTR 

T2bTR2*R2 

T3«TR3*R2 

T4«TR4*R2 

T5«TR5*R2 

EDT1«EXPJ-(T3-R2)/THETAR) 

EDT2«EXP<-(T3-T2)/THETAR} 

EDT3-EXP  <- (T5-T3) /THETAR) 

EDT4«EXP(.(T5-T4)/TMETAR) 

302  CALL  STPWB(T2|VMID.Zl.l6> 

CALL  STPW9(T3,STPW,Z1.16> 

FUFI*EDT1+(VMID*U*E0T2  /TR2*STPW*0«5*U/TR3>*DSU 

PRE-STPW 

T»T3 

N-5 

GO  TO  71 

305  CALL  STPWB,’T4,VMI0#Z1i16) 

CALL  STPWB (TStSTPW tZi f 16) 

FI»FI*EDT3* ( 0,5*PRE*U»EOT3/TR3*3,»VMlD*U*EDT4/TR4«STPW*U/TR5) *DSU 

PRE-STPW 

T-T5 

N«6 

GO  TO  71 
30B  N«13 

DTACT-2.*TACT*0T/3. 

310  IF ( (IPRES.LT.O) ,OR« (T.LT.TGU  GO  TO  MO 
320  DT-DT1 

DTACT«2.*DT*TACT/3. 

DST-DT/3. 

E0t-EXP(-DT/THETAR) 


CALCULATION  OF  MAIN  BOTTOM  REFLECTED  WAVE 

400  N»7 
410  T-T+DT 

CALL  STPWB(T.VMID.Zl.NP) 

EVERY  TENTH  STEP  RECALCULATE  STPW  WITH  HALF  THE  INTEGRATION 
POINTS  NP,  IF  ERROR  IS  LESS  THAN  .005  REDUCE  NP  BY  HALF. 
IF(MOD(IP»10).EO.O)  GO  TO  4lS 


BOTR602 
BOTR603 
BOTR604 
BOTR605 
BOTR606 
BOTR607 
ROTR608 
BOTR609 
BOTR610 
B0TR611 
BCTR612 
B0TR6I3 
BOTR6 14 
B0TR615 
B0TR616 
BOTR617 
B0TR618 
BOTR619 
BOTR620 
B0TR621 
B0TR622 
B0TR623 
B0TR624 
B0TR625 
R0TR626 
B0TR627 
B0TR628 
R0TR629 
BOTR630 
B0TR631 
B0TR632 
BCTR633 
B0TR634 
B0TR635 
B0TR636 
B07R637 
B0TR638 
B0TR639 
BOTR640 
B0TR641 
B0TR642 
B0TR643 
B0TR644 
B0TR645 
B0TR646 
B0TR647 
B0TR648 
B0TR649 
BOTR650 
BOTR651 
B0TR652 
BOTR653 
B0TR654 
B0TR655 
B0TR656 
B0TR657 
B0TR658 
B0TR659 
BOTR660 
B0TR661 
B0TR662 
B0TR663 
R0TR664 
B0TR665 
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IF (NP.EQ.4)  60  TO  415 

NP2aNP/2 

CALL  STPW0(T,VMIO2.Z1,NP2) 

ERR0R«A8S( (VMIl-VMI02) /VMJO) 

IF(ERR0R.LT. 0.005)  NP.NP2 
415  T*T*0T 

CALL  STPWB(T,STPW,Z1.NP) 

FlaFI*E0T**2*<  <PRE*£DT*4.*vMID)*E0T*STPW)*DST 

?:;c«stpw 

OO  TO  71 


CALCULATION  OF  0OTTOM  REFLECTION  FOR  TIMES  LESS  THAN  OR  EQUAL  TC 

the  arrival  time  of  the  peak  if  the  angle  of  incidence  is  less 
than  or  equal  to  the  critical  angle,  the  direct  nave  arrives 

0EFORE  OR  TOGETHER  NITH  THE  BOTTOM  REFLECTION. 

700  T-1.0 
Nag 

707  STPWaO, 

GO  TO  71 
706  N*U 
700  T»T*DT*2»0 

IF  <T.2.0*OT .LE.R2)  GO  TO  71 
710  T»R2 

STPW*CR/R2 

PRE*0. 

WRITE  (6 .550) 

NPag 

N»7 

00  TO  71 


CALCULATION  OF  OIRECT  WAVE  -PO-t  BOTTOM  REFLECTION  -PBOT-, 

SURFACE  REFLECTION  -PS-  AND  TOTAL  PRESSURE  -P-  . 

70  IF(T.LT.I.O)  GO  TO  72 

71  PD»PACT*EXP< (l.O-TJ/THETA) 

IF(T.LT.RS)  GO  TO  72 
PSa-PACT/Rl*EXP< (RS-T5/THETSR) 

72  P8OTlaP0OT 
FTHEtAcFI/THETAR 
PBOT=PaCT/PACTC* (STPN-FTHETA) 

73  TIME* TACT* (T-l oO) 

Pa  PO  ♦  PS  ♦  PBOT 

TEST  TO  INSURE  THAT  THE  ABSOLUTE  PRESSURE  CP.HYDROSTATIC)  .GE.  0. 

PaAMAXl (P.PH) 

IF(T.GT.TSTOP)  N=15 


CALCULATION  OF  IMPULSE*  FIMP  AND  ENERGY  FLUX*  EFLUX 
CALCULATION  OF  POSITIVE  IMPULSEa  POSIMP 

6004  XPaAMAXl  (P • 0 • > 

GO  TO  (6030.6007i6020.6030i604Q.6050i6007t6030.6007* 
I  6030.6010.6007,6007.6007.6060) .N 
6007  IF(IPRFS.LT.O)  GO  TO  6070 
PMIOaP 
XPM IO*XP 
GO  TO  6090 

6010  XXOT*TACT/3.*(T-TPRE) 

FIMPa3.*(P.PPRE)*XXDT  /2.*FIMP 


BOTR666 

B0TR667 

BOTR660 

B0TR669 

BOTR670 

B0TR671 

80TR672 

B0TR673 

B0TR674 

B0TR675 

B0TR676 

BOTR677 

B0TR678 

BOTR679 

BOTR6BO 

80TR6B1 

B0TR682 

BOTR683 

BOTR684 

B0TR685 

B0TR686 

B0TR687 

BOTR688 

B0TR689 

BOTR690 

B0TR691 

BOTR692 

B0TR693 

BOTR694 

B0TR695 

B0TR696 

B0TR697 

B0TR698 

B0TR699 

BOTR700 

BOTR701 

BOTR702 

BOTR703 

BOTR704 

B0TR705 

BOTR706 

3OTR707 

BOTR708 

BOTR709 

BOTR710 

B0TR7 1 1 

B0TR712 

B0TR713 

BOTR714 

BOTR715 

B0TR716 

BOTR717 

BOTR718 

BOTR719 

0OTR72O 

B0TR721 

B0TR722 

BOTR723 

B0TR724 

B0TR725 

B0TR726 

BOTR727 

80TR728 


POSIMP=3.* (XP.XPPRE) *XXDT/2 •♦POSlMP 


B0TR729 
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EFLUXwEFLUXMABS(P>*P*ABS<PPRF>*PPRE)*XXDT/2./RH0wAT/CWATER«3./ 

1  2(3066 

C  THE  CONVERSION;  FACTOR  2.3066  IS  NECESSARY  FCR  EFLUX  TO  HAVE 
C  UNITS  IN-PSI 
PPRE«P 

xppre*xp 
IPRES-1 
GO  TO  6092 

6020  OTACT«2.0*OSV«TACT 
PPRE«PPRE<*V/TR1 
XPpRE«XPPRE*V/TRl 
PMin«P*V/2,/TR3 
XPHlO«XP*V/2./TR3 
GO  TC  6090 
6030  PENO«0. 

xpend-o. 

GO  TO  6072 

6040  OTACT«2.0*DSU#TACT 
PPRE*0. 

XPPRE«0. 

PMin«P*U/2,/TR3 
XPMID«XP*U/2./TR3 
GO  TC  6090 
6050  peno«p*u/tr5 
XPEN0«XP*U/TR5 
GO  TO  6072 

6060  IF(IPRES#GT.O)  GO  TO  6010 
6070  PEND-P 

XPEKD-XP 

6072  FIMP-FIMP* <PPRE*4.«PMIC*PEnO) *DTACT 

P0SIHPaP0SIMP*(XPPRE*4.*XPMlD*XPEND)«0TACT 

EFLUX«EFLUX*(A0S(PPRE)#pPRE*4.*A8S(PMiD)»PMID*ABS(PEND)«PEND)» 

1  OT AC T/RHOW AT/C WATER/2 .3 066 
PPRE*P 
X°PRE«XP 

6090  IPf'FS«-l«IPRES 
6092  IF ( IPRES.GT . 0 )  TPRE«T 

WHEN  02R2.LT, 1,0  THE  BOTTOM  REFLECTION  IS  NOT  CALCULATED  AT 
THE  nIRECT  WAVE  ARRIVAL  TIME  T*l,0.  IN  ORDER  TO  PLOT  THE 
INSTANTANEOUS  RISE  OF  THE  DIRECT  SHOCK  AT  T«I,0t  THE  BOTTOM 

reflfction  is  obtained  by  linear  interpolation,  plot  points  ARE 

THEN  CALCULATED  FOR  THE  TCP  AND  BOTTOM  OF  THE  SHOCK  FRONT, 

I F { (IP.EO.l) .OR. (IP. GT. 1000) >  GO  TO  7002 
IF( (TIME. GT.O.) .ANO. (XX  UP-1) ,LT.O,))GO  TO  6095 
GO  T"  7002 
6095  X X ( I b ) >0 , 

IF (T.NE.R2)  P0OTO=P0OT1 ♦XXflP-lJ* (PBOT“PROTl >/(XX(IP-li"TlME*l,E6) 
IF (T.EC.R2)  P80TO*P90T1 
YY(Ip)»A*'AX1  ( proto »Ph; 

XX (Io*l)«0. 

yy(Io*i)*amaxhphotd*pact,ph) 

WRITp(6»559)  YY(IP*1> 

IP*Ib*2 

C  PRINT  ROUTINE 

C  FORMATS  are  LISTED  AT  THF  END  OF  THE  PROGRAM,  CARDS  ROTR907-1042 

7002  *F  (Z5.GT.0.0)  GO  TO  7003 

WRITE  (6»551 )  T»STPW,FTHETA,PD*TIMFtP.FIMP 
WRITE (6»59B)PR0T»PS,VMlDfPRE»EFLUX 
GO  TO  7004 

7003  WHITE (6»551 )  T»PBOTr EFLUX »PD* TIME ,P,F IMP 
GO  TO  7005 

C 

C  REDUCED  IMPULSE 

7004  RF IMP=FIMP/W13R/RACTU 


B0TR730 

B0TR731 

BOTR732 

R0TR733 

B0TR734 

B0TR735 

R0TR736 

B0TR737 

BOTR73B 

B0TR739 

BOTR740 

B0TR741 

B0TR742 

B0TR743 

B0TR744 

R0TR745 

B0TR746 

B0TR747 

B0TR748 

BOTR749 

BOTR750 

BOTR751 

B0TR752 

ROTR753 

B0TR754 

BOTR755 

BOTR756 

0OTR757 

BOTR7S8 

BOTR759 

BOTR760 

B0TR761 

BGTR762 

B0TR763 

B0TR764 

B0TR765 

B0TR766 

80TR767 

B0TR768 

80TR769 

BOTR770 

0OTR771 

BOTR772 

0OTR773 

B0TR774 

BOTR775 

0OTR776 

0OTR777 

80TR778 

0OTR779 

BOTR780 

BOTR78I 

BOTR782 

0OTR783 

0OTR7B4 

BOTR7B5 

BOTR786 

B0TR787 

BOTR788 

0OTR789 

HOTR790 

BOTR791 

B0TR792 

B0TR793 
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;  REDUCED  POSITIVE  IMPULSE 

RPOSIMaPOS IMP/W13R/RACTU 
REDUCED  ENERGY  FLUX 

SriteigUsm^esid^fihp.reflux.posimp.rposim 
7005  IF(IP.GT.IOOO)  GO  TO  999 

XX<IP)cTlME*I.0E6 

YY { IP) =P 

:  BOTTO^REFLECTION  TIME  AND  PRESSURE  STORED  IN  OX  AND  GY  FOR 

C  PTV  calculation. 

IF (T.LT.D2R2)  GO  TO  7001 
IPTV=IPTV*1 

OX { I PTV) =TACT* (T-D2R2 ) 

QY (IPTV)*AMAX1 (PBOT.PH) 

GO  TO  (H0*1 59 *210* 30 0*3 05 *308 *410*706 *201*159 *159 *114* 

1  310*708*999) »N 


7001 

: 

999 


IPMAX*IP“1 

calcomp  plot  tape  generated  if  Z3  *  o* 

CALL3PL0T1(Xx!yyIiPMAX.X1.X2.ADATE.THE*WCH.CB0T.P0ISR,Z1) 

;  997  CALCULATION* OF  INPUTS  FOR  SUBROUTINE  PTV. 

IF (RADIUS.LE.Ot )  GO  TO  4 
1997  TIMER2»TACT*(R2“D2R2) 

XT3«TIMER2-2.»DT*TACT 
XT4»0.8*TIMER2 
XT5bQX ( IPTV) 

IPTV»IPTV*1 
QX  < IPTV) “1 »0E20 
QY(IPTV)*0. 

COSA*COSTH*COS(A)*SINTH*SIN(A)  cl„rTTV  IDTU, 

r  CALCULATION  OF  PEAK  TRANSLATIONAL  VELOCITY  (PTV). 

CALL  PTV (TIMER2.XT3*XT4.XT5.RADIUS,30..APRINT.COSA,RHOWAT. 

1  CWATER.TIME1.PTV1.PTV2) 

IF (APR I NT .GT.O. )  WRITE16.570)  TIMEl.PTVl.PTV2 
GO  TO  4 


PLANE  WAVE  APPROXIMATION  USING  EQUATIONS  OF  ARONS  AND  YENNIE 


998  IF (CSHEAR.GT .0. )  GO  TO  1005 
WRITE (6.590) 

GO  TO  1007 

1005  FORMATS* ARE* LISTED  AT  THE  END  OF  THE  PROGRAM.  CARDS  BOTR907 
SELECTION  OF  TIME  STEP 

1007  WRITE (N.591 ) 

IF  (DTR2.GE.R2)  GO  TO  1020 

l  CALCULATION  OF  TIME  STEP  FOR  SUPERCRITICAL  REFLECTION 

1010  Ms  (R2-D2R2) *STEPS/THETA/2.*1.0 
M«2*M-1 

IF  (M.LE .4)  GO  TO  1020 
1012  DTs«R2-D2R2) /FLOAT (M) 


-1042 


B0TR794 
B0TR795 
B0TR796 
B0TR797 
B0TR798 
B0TR799 
80TR800 
BOTR801 
BOTR802 
BOTR802A 
BOTR802B 
B0TR802C 
B0TR802D 
BOTR802E 
B0TR8C2F 
BOTR803 
BOTR804 
BOTR805 
BOTR806 
BOTR807 
BOTR808 
BOTR809 
BOTR810 
BOTR811 
BOTR812 
BOTR813 
B0TR813A 
8QTR813B 
B0TR813C 
B0TR8130 
80TR813E 
B0TR813F 
B0TR813G 
B0TR813H 
B0TR813I 
B0TR813J 
B0TR813K 
B0TR813L 
B0TR813M 
B0TR813N 
BOTR814 
B0TR815 
»«  B0TR816 
BOTR817 
BOTRP18 
BOTR819 
BOTR820 
BOTR821 
BOTR822 
B0TR823 
B0TR824 
BOTR825 
B0TR826 
BOTR827 
B0TR828 
B0TRP29 
BOTRP30 
B0TR831 
B0TR832 
B0TR833 
BOTR834 
B0TR835 
B0TR836 
B0TRP37 
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IF (D2R2-1.) 1014 ,1015. 1015 

1014  T»02R2 

GO  TO  1700 

1015  T«l. 

GO  TO  1700 

CALCULATION  OF  TIME  STEP  FOR  SUBCRITICAL  REFLECTION 
1020  DTcTHETA/STEPS 
T«l. 

1700  IF(T.GE.l.O)  GO  TO  1731 

1730  SW=0. 

GO  TO  1732 

INCIDENT  (DIRECT)  WAVE  RESPONSE 

1 731  SW«EXP(-(T-1.0)/THETA) 

1?32  XE1=0. 

XEI=0. 

IF«E2.NE.O.)  GO  TO  1738 
1733  IF(T.GE.«2)  GO  TO  1736 
1*35  PRFL=0. 

GO  TO  1745 

PRESSURE  RESPONSE  FOR  SUBCRITICAL  REFLECTION 
1736  PRFL«CR*EXP<-<T-R2>/THETAR) 

GO  TO  1745 

1738  TBTH*(T-R2)/THETAR 

IF(TBTH)  1741,1742.1743 

1741  CALL  EXEl ( TBTHt XE1 ) 

PRESSURE  RESPONSE  FOR  PRECURSOR  OF  SUPERCRITICAL  REFLECTION 
PRFL*CR  'SIN(E2)*XE1/PI 
GO  TO  1745 

PRESSURE  RESPONSE  AT  SINGULARITY 

1742  PRFL»(1.E*30)*SIGN(1..EE) 

GO  TO  1745 

1743  CALL  ExEI(TBTH.XEI) 

PRESSURE  RESPONSE  FOR  MAIN  SUPERCRITICAL  BOTTOM  REFLECTION 
PRFL*CR*(EXP(-TBTH)*C0S(E2)-XEI*SIM<E2)/PI) 

1745  PRFL»PRFL/PACTC/R2 

pbot«pact«prfl 

1 7 1 0  TIME=TACT«(T-1.) 
po»pact«sw 

IFd.GF.RS)  PS*-PACT/R1*EXP<  (RS-T)/THETSR) 
total  PRESSURE  (NEGATIVE  VALUE  LIMITED  TO  PH) 

P«AMAX1 (PD+PBOT+PS.PH) 

OUTPUT  ROUTINE 

WHITt r 6.551 )T,PBOT.PD. PS. TIME. P 

IF(IP.GT.IOOO)  GO  TO  4000 
7100  XX<IP)=TlME«1.0E6 
YY  < IP) sP 
IP«IP+1 

BOTTOM  REFLECTION  TIME  AND  PRESSURE  STORED  IN  QX  AND  QY  FOR 
PTV  CALCULATION, 

IF(T.LT.D2R2)  GO  TO  7102 
IPTV=IPTV+1 

QX  < IPTV) *TACT* (T-D2R2) 

QY  < IPTV) *AMAX1 (PBOT.PH) 

7102  TDT»T 
T»T*DT 

IF (ABS (TDT-R2) ,LT.1,5*DT)  T=TDT*0,2*DT 
I F  < (TDT.LT.R2) .AND. (T.GT.R2) )  T=R2 


BOTKP38 

BOTR839 

BOTR840 

B0TR841 

8OTR042 

BOTR843 

80TR844 

B0TR84E 

B0TR846 

BOTR847 

90TR848 

BCTR849 

EOTR850 

BOTR851 

30TR852 

BOTR853 

BOTR854 

30TR855 

80TR856 

80TR857 

B0TR858 

30TR659 

BOTRS60 

B0TR861 

BOTR862 

BOTR863 

B0TR864 

B0TR86S 

B0TR866 

BOTR067 

B0TR868 

B0TR869 

B0TR870 

B0TR671 

B0TR872 

BOTR873 

BOTR874 

BOTR875 

BOTRP76 

BOTR877 

BOTR878 

BOTR879 

B0TR880 

BOTR881 

BOTR082 

B0TR883 

B0TR884 

BOTR085 

BOTR886 

BOTR887 

BOTR888 

BOTR889 

BOTR890 

B0TRS91 

B0TR891A 

B0TR8918 

B0TRB91C 

BOTR091D 

B0TP891E 

B0TR891F 

80TR892 

BOTR893 

B0TR894 

BOTRP95 


A-15 


o  o  o 


r*-r  '•  fS.m  lV" 


jPWV’P’  I 


NOLTP  71 -110 


IF( <TDT,LT.1.0> .AND. <T.GT.1,0) )  T=1.0  80TR896 

IF(T.LT.TSTOP)  GO  TO  1700  80TR897 

4000  IPMAX»IP-1  80TR898 

C  CALCOMP  PLOT  TAPE  IS  GENERATED  IF  73  >  0,  B0TRG99 

IF(Z3,NE*0.)  GO  TO  4001  80TR900 

CALL  PLOTl <XX,YY.IPMAX»X1.X2.ADATE.THE»WCH»C80T»P0ISR,Z1)  80TR901 

4001  ICASE*ICASE*1  BOTR902 

IF (RADIUS, GT.0.)  GO  TO  1997  B0TR902A 

GO  TO  4  6OTR903 

80TR904 

***********************************<**«************************<****  80TR905 
„  80TR906 

500  FORMAT  (10X.26HSLOPE  OF  BOTTOM  IN  DEGREES. 54X.6HSL0PE*.E16, 8  )  BOTR907 

501  FORMAT (EX, 12HREDUCED  TIME.23X.11HENERGY  FLUX.26X.4HTIMp.llX.  8OTR908 

1BHPRESSURE.10X.7H  IMPULSE/ 1  OX. 1HT.16X.4HPB0T.31X.2HPD.12X,  BOTR909 

27HSEC0NDS.14X.3HP Si  /  )  BOTR910 

502  FORMAT (10X.39HREDUCED  SLANT  DISTANCE  (RACTU/WCH*«l/3 ) .41X.6HREDR*  80TR911 

1.1E15.8)  80TR912 

503  FORMAT (10X.60HPRINT  OUT  CONTROL  PARAMETER  (Z5.GT,0.  FOR  SHORTER  PR  BOTR913 

1INT  OUT) »20X.3HZ5*.E16,6)  B0TR914 


504  FORMAT (10X.40HWEIGHT  OF  EXPLOSIVE  CHARGE  IN  LB  (OR  KT) .40X.5HW0H*  B0TR915 

1  E15.B  )  80TR916 

505  FORMAT (10X.36HVELOCITY  OF  SOUND  IN  WATER  IN  FT/SEC. 44X.8HCWATER*  .  80TR917 

1E15.8  )  80TR91B 

506  FORMAT (10X.37HVELOCITY  OF  SOUND  IN  BOTTOM  IN  FT/SEC. 43X.6HCB0T*  .  80TR919 

1E15.8  )  80TR920 

507  FORMAT (10X.25HDENSITY  OF  WATER  IN  GM/CC.55X.BHRH0WAT*  ,E15.B  )  80TR921 

508  FORMAT (10X.26HDENSITY  OF  BOTTOM  IN  GM/CC.54X.BHRH0B0T*  .ElS.B  )  BOTR922 

509  FORMAT(10X.51HDURATION  AFTER  DIRECT  ARRIVAL  IN  MULTIPLES  OF  THETA,  BOTR923 

129X,7HDURAT=  1E15.8  )  B0TR924 

510  FORMAT ( 1H1 ,52X . 17HBOTTOM  REFLECTION. 10X.4HDATE.2X, 1A10  )  B0TR925 

511  FORMAT (10X.20HDEPTH  OF  WATER  IN  FT.60X .6HBIGH*  1E15.B  )  80TR926 

512  FORMAT(10X.20HDEPTH  OF  GAUGE  IN  FT.60X.6HDGAU*  1E15.8  )  BOTR927 

513  FORMAT (10X.24HDEPTH  OF  EXPLOSION  IN  FT.56X.3HD*  1E15.8  )  80TR928 

514  FORMAT (10X.50HHORIZONTAL  DISTANCE  BETWEEN  CHARGE  AND  GAUGE  IN  FT. 3  B0TR929 

10X.8HSMALLR*  1E15.8  )  BOTR930 

515  FORMAT (10X.41HCOEFFICIENT  OF  SW  PRESSURE  FORMULA  IN  PSI.39X,  BOTR931 

1  BHPRECOE*  .E15.8  )  BOTR932 

516  FORMATUOX.31HEXPONENT  OF  SW  PRESSURE  FORMULA, 49X.BHPREEXP*  .1E15,  80TR933 

18  )  80TR934 

517  FORMAT (10X.50HCOEFFICIENT  OF  SW  TIME  CONSTANT  FORMULA  IN  SECONDS.  80TR935 

1  30X.BHTHECOE*  .E15.8  )  BOTR936 

518  FORMAT (10X.36HEXPONENT  OF  SW  TIME  CONSTANT  FORMULA. 44X.8HTHEEXP*  *  BOTR937 

11E15.8  )  B0TR938 

519  FORMAT (10X.31HNUM8ER  OF  SUBDIVISIONS  OF  THETA  49X.7HSTEPS*  .1E15.B  BOTR939 

1  )  BOTR940 

520  FORMAT (1H0.47X.25HCHARACTERISTIC  MAGNITUDES  /  )  B0TR941 

521  FORMAT (10X.33HANGLE  OF  INCIDENT  WAVE  IN  DEGREES. 47X.5HTHE*  .1E15.8  BOTR942 

1  )  B0TR943 

522  FORMAT (10X.45HCRITICAL  ANGLE  OF  COMPRESSION  WAVE  IN  DEGREES. 35X.7H  B0TR944 

1ALPHA*  1E15.8  )  B0TR945 

523  FORMAT (10X.33HREDUCED  TIME  OF  PRECURSOR  ARRIVAL, 47X.6HD2R2*  .IE15.  BOTR946 


IB  )  807R947 

524  FORMAT (1  OX. 35HREDUCED  TIME  OF  GROUND  WAVE  ARRIVAL. 45X.4HR2*  .1E15,  B0TR948 

IB/  )  B0TR949 

525  FORMAT (10X.68HSLANT  DISTANCE  BETWEEN  CHARGE  AND  GAUGE-CHARACTERIST  80TR950 

1IC  LENGTH  IN  FT.12X.7HRACTU  *  .1E15.8  )  80TR951 

526  FORMAT (10X44HCHARACTERISTIC  TIME*RACTU/CWATER  IN  SECONDS, 36X.6HTAC  BOTR952 

IT*  , 1E15.8  )  BOTR953 

527  FORMAT (10X.58HCHARACTERISTIC  PRESSURE*FREE  WATER  SW  PEAK  PRESSURE  B0TR954 

1 IN  PSI  22X.6HPACT*  .1E15.8)  B0TR955 

528  FORMAT (10X.38HREDUCED  TIME  CONSTANT  OF  INCIDENT  WAVE.42X.7HTHETA*  B0TR956 

1.1E15.8  )  BOTR957 

529  FORMAT  (1HD  80TR958 
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-  -4X.17HSTF.PWAVE  RESPONSE. 3X.1?HC0NV0LUTI 

S3°XON-':lsi.*5T1I«”ll*%HP^SSURE.JiOK.,H,«P»L|E/aMpo>lEX.THSECON1JS  BOmM 

z  f[Jx/6111i0„seco™  bo».  ROB.JX.  7M 

,,2%m!iTa”rE°  „f  PEAK  OF  BOTTOM  REELECTED  »**E.  “’"oSI 

533  FORMAT (  IOXxAShREDOCEO  TIME  BOTRR6S 

>  OF  PRECURSOR  /  I  n/  „,5„SMALLH.UX.B"0EZE  BOTRRn 

535  fSS«T!Ux.aWM«TJ«S,fMT>«c'5TSA8,.SHSIBTH/  >  OOTRRT2 

»,SiVKT,°  8ET“EEN  ,NC,nENT  ZZ  KSSi 

5371TICALTanOLE  IS  TMOXAE/t  betrEEN  1NCI0ENT  AMO  CRITICAL  AN  •  -  BO,„,76 
538 ME  « START  «  R--O-.A0X.AH  TMET= 

-X.—  — «*  ream,,  gg 

sp»sr  si 

■  - 1 
s 

FORMATtlox.5lHBOTTOM  REFLECTED  w*vt  B0TR994 

™‘riS5«SSSTW*^  •»■'  5”-5"INP . TIPE  0F  SROUno..»e  M 

a CHM0E0  50  TH" ARR,V  _  ssmseore  amo  Sks 

—  "UI RE5ULTS  SS 

a  MAY  BE  MEANINGLESS  //  >  BOT1003 

in  rsss  sri  .  -  -  ski:k 

55<j  FORMAT  (IHObIOX.P  n0n-r1G1D  BOTTOM  )  B2iinn7 

ii?  ;s“;ss,»:i,.Sto-  »a»e  >  Un!.. 

{**  FfSSTT!K,.:Sa,.!S  KnO«  i'.TH  SET.  SHEAR  .«*  ,  |gilS?0 

c  565  FORMAT  |4BX,?*HRL“EE«AVEi«PPMX>"”|Hn  ,  B0T1MJ 

566  E°HMftT 'll!’ iIhROSENBAUM  METHOD  >  5RX ,TMRADlOS=  »E15*S  >  HOT  in 1  IE 

Si 

1E15.B  .  10H,.... . .  GEOMETRY  CHANGE  ROT,. IT 

574  FORMAT (5X. 10n* 

<  nuTKir.  nnTTOM  ) 
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579  FURMAT (10X.55HPEDUCED  ARRIVAL  TIME  OF  CRITICALLY  REFRACTED  SHEAR  W 
1 AVE.25X.8HSHD2R23  EI5.8) 

584  FURMAT (10X,52HSCALIN6  PARAMETER  FOR  Y-AXIS  (PSl  PER  INCH  OF  GRAPH) 
1  2BX  *  4HX 1  a  ,E15.8  ) 

535  FORMAT  ( 1  OX » 61HSCALING  PARAMETER  FOR  X-AXIS  (MICROSECONDS  PER  INCH 
10F  GRAPH)  19X  »4HX2»  »E15.8  ) 

586  FORMAT (10X.29HPARAMETER  THAT  SELECTS  TNE0RY«51X»4HZ1»  .E15.8  ) 

587  FORMAT  (I0X.43HARRIVAL  TIME  OF  GROUND  WAVE  IN  MICROSECONDS  37X.4HZ2 

1=  »E15 ,8  ) 

588  FORMAT (10X.55HPLOT  CONTROL  PARAMETER  (Z3  ■  0.  MEANS  PLOTS  ARE  WANT 
1ED)  1 25X • 4HZ3a  »E15.8) 

590  FORMAT  (43X.39HAR0NS-YENNIE  APPROACH  NON-RIGID  BOTTOMS  /  ) 

591  FORMAT  (  5X» I2HREDUCED  TIME.3X.17HBOTTOM  REFLECTI0N»4X,9HSH0CKWAVE» 

1  4X » 18HSURFACE  REFLECT  ION .5 X »4HTI ME »8X , 14HT0TAL  PRESSURE/lOXt 

2  IHT,15X»4HPB0T»14X.2HP0*15X»2HPS.)2X*7HSEC0NDS*9X* 

3  7HP  (PSD  /  ) 

592  FORMAT (10Xt22HREFLECTION  COEFFICIENT. 58X,4HCRa  1E15.B  ) 

593  FORMAT (10X.30HANGLE  OF  PHASESHlFT  IN  DEGREES  »50X»*HEEi  1E15.8) 

594  FORMAT (10X»39HANGLE  OF  SHEARWAVE  IN  BOTTOM  IN  DEGREES »41X »6HANGA* 

I  lE15.fi  ) 

596  FORMAT (37X.47HAR0NS-YENNIE  APPROACH  EXTENDED  TO  RIGID  BOTTOMS  /  ) 

597  FORMAT (10X.43HANGLE  OF  PRESSURE  WAVE  IN  BOTTOM  IN  DEGRfES.37X.0FTH 
1EONE*  •  1E15.8  ) 

598  F0RMAT(18Xt4El7.7,17XtE!7.7) 

599  FORMAT (10X.42HINPUT  INCONSISTENT.  COMPUTATION  SUPPRESSED/  ) 


END 


•  ••••iJURROUTINE  STONL**** 


CALCULATION  OF  PROPAGATION  VELOCITY  of  STONLEr  WAVE, 

SUBROUTINE  stonl 

common  B#C0SALtC0STHtR2,SlNBE,SlNTHtCwATER»CB0T.CSHEAR.CST0N.PESlD 
COMMON  C2,CBOT2#CSHR2,C8SH,C2SHR2,C4C0,SINTH2 

IF(CSHEAR.lE,0.)  GO  TO  4l6 
4  Cl2aCwATER«*2 
C32aCBCT**2 
C4?*CSHEAR##2 

iteration  process 
IF(CWATER,GT . CSHEAR)  go  to  2 

1  Y2aCWATFR##2 
GO  TO  3 

2  Y?*CSHEAR«#2 

3  CKsy2/1000. 

FYaSCRT (C12-Y?)*(CB0T« (y2-2.*CA2)**2-4,*CSHEAR*C42*SQRT<  <C32-Y2)»( 
1C42-Y2 ) ))*B«CWATER*Y2**?*S0RT(C32-Y?) 

Y2aY2-CK 

400  00  410  IR«1 .999 

FS  STORED 
FSaFY 


STOKLEY  wave  VELOCITY 

FYaSQRT (C12-72) *(CBOT* (Y2-?,*CA2)**?-4,<»CSHEAP#C42*SQRT ( (C32-Y2*  * ( 
1C42-Y2) ) ) ♦B«CwATER*Y2*«2*S0RT (C32-Y2) 


BOT 10 18 
BUT  10 19 
80T1020 
BOT1021 
BOT1022 
BOT1023 
BOT1024 
BOT1025 
BOT1026 
BOT1027 
BOT 1028 
BOT1029 
BOT1030 
BOT1031 
BOT1032 
BOT1033 
BOT1034 
BOT1035 
BOT1036 
BOT1037 
30T1038 
BOT1039 
BOTIOAO 
BOT1041 
8OT1042 
BOT1043 
BOT10A4 


STONOOl 

STON002 

STON003 

STON004 

STON005 

STON006 

STON007 

STON008 

STON009 

STONOlO 

STONOU 

STON012 

STONOU 

STONOU 

STONOU 

STONOU 

STON017 

STON018 

STONOU 

STON020 

STON021 

STON022 

STON023 

STON024 

STON025 

STON026 

STON027 

STON028 

STON029 

STON030 

STON031 

STON032 

STON033 

STON034 
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STON035 

IF(FY)412«4l5,408  STON036 

408  Y2»Y2-CK  STON03? 

410  CONTINUE  STON038 

WRITE  <6*401 )  CWATER, CBOT. CsHEaR,B,Y2,FS.FY  STON039 

STOP  STON040 

ST0N041 

FV  IS  NEGATIVE  STON042 

STON043 

FALSE  POSITION  OR  SECANT  METHOD  ITERATION  STON04* 

412  YS«Y2*CK  STON04S 

DO  450  1*1 »50  STON046 

YSS*Y2  STON04? 

IF(ABS((YS-V2)/YS).LT.1.0E-7)  GO  TO  4i5  STON048 

440  Y2*YS*FS*(Y2-YS)/(FS-FY)  STON049 

PS*FY  STON050 

YS.YSS  STON051 

STON052 

FYuSQRT (C12-Y2) * (CBOT* ( Y2-2«*C42) ••2-*,#CSHEAR#C42*S0aT { (C32-Y2) * (  STON053 
1C42-Y2) ) ) ♦B*CWATER*Y2**2*SQRT <C32-y2>  STON054 


STON055 

450  CONTINUE  STON056 

WRITE (6*402)  CWATER* CBOT tCSHEARtB* y2 *YS*FS»FY  STON057 

STOP  ST0NC5S 

STON059 

RESULT  STON060 

415  CST0N«SQRT(Y2)  STON061 

RETURN  STON062 

416  CST0N*-0.  STON063 

RETURN  STON064 

STON065 

401  FORMAT (20X.  42HFIRST  ITERATION  FOR  CSTON  DID  NOT  CONVERGE//  STON066 

1  30H  CWATER. CBOT. CSHEAR,B,Y2.FS*FY  //  1P7E16.6)  STON067 

402  FORMAT (20X.  43HSEC0ND  ITERATION  FOR  CSTON  DID  NOT  CONVERGE//  STON068 

1  33H  CWATER. CBOT. CSHEAR.B . Y2. YS.FS.FY  //  1P8E14.6)  STON069 

STON070 

END  STON071 


»**»*SUBROUTINE  stpwa**** 


STPA001 


STPA002 

STPA003 


PRECURSOR  CALCULATION  USING  CAGNIArD  METHOD.  STPAOO* 

STPA005 

STPA006 

SUBROUTINE  STPWA(T.STPW.CONTR.K)  STPA007 

DIMENSION  P (30)  STPA008 

COMMON  B.COSAL. C0JTH.R2.SINBE.SINTH. CWATER sCBOT.CSHEAR. CSTON. RESID  STPA009 
COMMON  C2.CB0T2.CSHR2.CBSH.C2SHR2.C4CB.SINTH2  STPAOlO 

EXTERNAL  ONE. SEVEN  STPA011 

DATA  AA.BB.SQ2/-1. 57079633, 1.57079633, 1.41421356/  STPA012 

TR*T/R2  STPA013 

V=SQRT(1.-TR**2)  STPA014 

5  IF (CONTR.EQ.3. )  GO  TO  100  STPA015 

STPA016 

STPA017 

CALCULATION  OF  THE  PRECURSOR  USING  CAGNIARD-'VOSENBAUM  INTEGRALS,  STPAOiB 

STPA019 

P  (9 ) so  ,  STPA020 

XM»COSTH«TR*SINTH*V  STPA021 

P  ( 1 ) sCOSAL-XM  STPA022 

P(2)»4.*V*SINTH/P(1)  STPA023 

P(5)=C0SAL*XM  STPA024 
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FACT0RrSQ2»B/ <R2»BB ) 


S7PW«FACT0R*FGI <AA»BB»K»ONE,P) 
RETURN 


CALCULATION  OF  THE  PRECURSOR  USING  COMPLEX  ARITHMETIC  METHOD 

100  P(1)>0. 

P (2) sTR 

P  (3 ) sCWATER/CBOT 
P  (4) «SINTH*TR-COSTH*V 
FACTOR«R2*BB 

ANS2*SEVEN(P(4) »P) 

STPWa(ANS2*FGI (P (3) »P<4) »K,SEVEN*P) )/FACTOR 

RETURN 

END 


STPA025 

STPA026 

STPA027 

STPA028 

STPA029 

STPA030 

'STPA031 

STPA032 

STPA033 

STPA034 

STPA035 

STPA036 

STPA037 

STPA038 

STPA039 

STPA040 

STPA041 

STPA042 

STPA043 

STPA044 

STPA045 

STPA046 


•SUBROUTINE  STPMB*»»« 


CALCULATION  OF  "'IN  BOTTOM  REFLECTION 


SUBROUTINE  STPWB (T «STPW«CONTR*K) 

DIMENSION  POO) 

COMMON  B»COSAL«COSTH#R2,SlNBEtSINTHtCMATER.CBOT.CSHEAR,CSTON.RESID 

COMMON  C2.CB0T2.CSHR2«CBSH,C2SHR2»C4CB.SINTH2 

EXTERNAL  TWOtONEl .ONEtSEVEN 

DATA  BP., SQ2/1. 57079633, 1.41421356/ 

TR-T/R2 

5  IF (C0NTR,E0.3. )  GO  TO  100 

CALCULATION  OF  THE  MAIN  REFLECTION  USING  CAGNIARD-ROSENBAUM 
INTEGRALS. 

P<9).  ,. 

MAGNITUDES  K  AND  L 
XK«COSTH»TR 
XL»SINTH*»2» (TR**2-1 * J 

P (7) aXK 
P(8)«XL 

MAGNITUDES  D»E»  and  F 
P  ( 1 1 ) »TR**2* ( 1 .-2,»SINTH**2) *SINTH»*2 
P (12) «4 ,»S INTH**2*C0STH**2»TR*»2» ( TR**2-1 # ) 

P  ( 13) »TR»*2-SINTH**2 
FACTOR*B/BB/R2 
IF(CSHEAR.GT.O.)  GO  TO  12 

NON-RIGID  bottoms 


10  TERM1»<1.-B)/(1.*B)/R2 

resid*o. 


STPB001 

STPB002 

STPB003 

STPB004 

STPB005 

STPB006 

STPB007 

STPB008 

STPB009 

STPB010 

STPBOU 

STPB012 

STPB013 

STPBOU 

STPB015 

STPBOU 

STPB017 

STPBOIS 

STPBOU 

STPB020 

STPB021 

STPB022 

STPB023 

STPB024 

STPB025 

STPB026 

STPB027 

STPB028 

STPB029 

STPB030 

STPB031 

STPB032 

STPB033 

STPB034 

STPB035 

STPB036 

STPB037 

STPB038 

STPB039 
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IF(CBOT.OT.CWATER)  00  TO  11 
SLOW  NON-RIGID  BOTTOMS 
SI6M«SQRT ( (CWATER/CBOT) **2-1 . ) 
STPWsTERMl-S02*FACT0R*F0I (0. .SIGM,K»TWOtP> 
RETURN 

FAST  NON-RIGID  BOTTOMS 
11  STPW«TERM1*FACTOR*FGI(O..COSAL»K»OnE»P) 
RETURN 


RIGID  BOTTOMS 


STONELEY  POLE  RESIDUE 

12  TERM«1,/R2 

CWS2«(CWATER/CSHEAR)**2 

SK«CWATER/CSTON 

RK2«SK#*2 

XG1«SQRT (ABS (SK2-1 • ) ) 

XG3«S0RT (ABS (SK2- (CWATER/CBOT) «*2> ) 

XG4=SQRT (ABS (SK2-CWS2) ) 

XSA«R2**2* (TR**2-SK2*C0STH**2) 

XSF« (2.*R2**2#TR*C0STH*XG1 ) ##2 

XNUM«XG1*( (CWS2/2.-SK2)**2-SK2*XG3*XG4)-B*X33«CWS2*»2/4. 
XDEN«((CWS2/2.-SK2)**2-SK2*XG3*XG4)/XGl-XGl*(2.*CWS2-4.*SK2*2.»XG3 
1*XG4^SK2* (XG4/XG3+XG3/XG4) ) *B*CWS2**2/4./XG3 

RESID*  -SQ2  *SQRT(ABS( (SQRT(ABS(XSA**2*XSF) )-XSA)/(XSA*»2*XSF) ) 
l)*XNUM/XDEN/XGl 

TERM1«1./R2*RESID 

IF(CSHEAR.GT.CWATER)  GO  TO  50 


SLOW  SHEAR  WAVE 
30  SIG2=SQRT (CWS2-1 • 1 

STPW»TERM1 ♦FACTOR* (FGI ( 0 . »C0SAL .K.ONE  » P ) -CWS2«*2*SQ2/4.» 
1  FGI (0 ( «SIG2«K tONEl »P) ) 

RETURN 

FAST  SHEAR  WAVE 

50  STPW«TERMl+FACTOR*FGI (0.»COSAL»K»ONEtP) 

99  RETURN 


CALCULATION  OF  THE  MAIN  REFLECTION  USING  COMPLEX  ARITHMETIC  METHOD 

100  P(l)oC0STH*SQRT(TR**2-l.) 

P (2) «TR 
P(3)«0, 

P(4)»SINTH*TR 

FACT0R«R2*BB 

ANS2*SEVEN (P (4) tP) 

STPW» ( ANS2*FGI (P (3) »P(4) tK,SEVEN»P) ) /FACTOR 


RETURN 

END 


STPB040 
STPB041 
STPB042 
STPB043 
STPB044 
STPB045 
STPB046 
STPB047 
STPB048 
STPB049 
STPB050 
STPB051 
STPB052 
STPB053 
STPB054 
STPB055 
STPB056 
STPB057 
STPB058 
STPB059 
STPB060 
STPB061 
STPB062 
STPB063 
STPB064 
STPB065 
STPB066 
5TPB067 
STPB068 
STPB069 
STPB070 
STPB071 
STPB072 
STPB073 
STPB074 
STPB075 
STPB076 
STPB077 
STPB078 
STPB079 
STPB080 
STPB0B1 
STPB082 
STPB0B3 
STPB084 
STPB085 
STPB086 
STPB087 
STPB088 
STPB089 
STPB090 
STPR091 
STPB092 
STPB093 
STPB094 
STPB095 
STPB096 
STPBC97 
STPB098 
STPB099 
STPB1 00 
STPB101 
STPB102 
STPB1C3 
STPB1C4 
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function  one««««  une 

ONE 

ONE 

integrand  of  the  cagniaro-rosensaum  integral  for  all  fast  bottoms.  one 

EXCEPT  FOR  PART  OF  THE  MAIN  REFLECTION  OF  A  BOTTOM  WITH  SLOW  SHEAR  ONE 

ONE 

ONE 


function  onecx.pj  one 

DIMENSION  P (30 )  ONE 

COMMON  B.C0SAL»C0STH.R2.SIN8E,SINTH.CWATER.C80T.CSHEAR»CST0N»RESI0  ONE 
COMMON  C2«CB0T2iCSHR2,CBSH»C2sHR2«C4CBtSINTH2  ONE 

ONE 

P  ARRAY  CALCULATED  IN  SUBROUTINES  STPwA  AND  STPW0  ONE 

TEST  FOR  PRECURSOR  PHASE  ONE 

IF(P(9)«GT.0«)  GO  TO  2  ONE 

ONE 

PRECURSOR  PHASE  ONE 

ONE 

SINX«SIN(X)  ONE 

W-(SINX«P(l)*P(S>>/2«  ONE 

XC2jC0SAL**2«W*«2  ONE 

IF(XC2.LT,0.)  XC2"0.  ONE 

FX*(1««SINX) •SORT ( ( (COSAL«W ) «P ( I ) ) / (1 , *S INX^P  ( 2 ) ) )  ONE 

GO  TO  3  ONE 

ONE 

MAIN  REFLECTED  WAVE  ONE 

2  Wax  ONE 

XC2aCOSAL*«2-W*«2  ONE 

IF(XC2.LT.O.)  XC2aO.  ONE 

RTlaSORT (XC2/ (P (8)  ♦ (W-P (7) ) ••2)1  ONE 

RT2aSQRT  (XC2/  (P  (8)  ♦  (W.P  (7) )  ••«  >  ONE 

FXaRTl-RT2  ONE 

ONE 

RELATIONS  FOR  PRECURSOR  AND  MAIN  WAVE  ONE 

ONE 

3  IF(CSHEAR.EO.O.)  GO  TO  110  ONE 

ONE 

RIGIO  BOTTOMS  ONE 

CSWaCSHEAR/CWATER  ONE 

CSW2"CSW**2  ONE 

PRC9«CSN2«(W««2-1,)*1,  ONE 

ONE 

XAaW*(l*"2«*CSW2*(l«"Waa2l )**2  ONE 

XBe4.«W«CSN2''CSW«(l.-W««2)«SQRT<XC2«A8S(PRCS) )  ONE 

ONE 

IF(FRCS.GE.O.)  GO  TO  22  ONE 

ONEaFX* (XA-XB) / ( ( X A»X8 ) ••2*B*B*XC2  )  ONE 

RETURN  ONE 

ONE 

22  XCeB«SQRT(XC2)  ONE 

ONEaFX*XA/(XA««2«(X8«XC)*«2)  ONE 

RETURN  ONE 

ONE 

NON-RIOIO  BOTTOMS  ONE 

110  ONEaFX*W/ (W»*2*B*B*XC2  )  ONE 

RETURN  ONE 

C  ONE 

END  ONE 


001 

002 

003 

004 

005 

006 

007 

008 

009 

010 

Oil 

012 

013 

014 

015 

016 

017 

018 

019 

020 

021 

022 

023 

024 

025 

026 

027 

028 

029 

030 

031 

032 

033 

034 

035 

036 

037 

038 

039 

040 

041 

042 

043 

044 

045 

04ft 

0^7 

048 

049 

050 

051 

052 

053 

054 

055 

056 

057 

058 
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C****«FUNCTICN  0NE1**** 

C 

C 

C  INTEGRAND  OF  THE  SECOND  CAGNI ARD-ROSEnBAUM  INTEGRAL  OCCURRING  FOR 

C  THE  MAIN  reflection  OF  A  bottom  WITH  A  SLOW  SHEAR  WAVE. 

C 

c 

FUNCTION  ONEl(X.P) 
dimension  POO) 

COMMON  8.C0SAL.C0STH.R2,SInBE,SINTH.CwATER.CB0T.CSHEAR.CST0N.RESID 
COMMON  C2.CB0T2,CSHR2,CBSH,C2sHR2,C*CB.SINTH2 
C 

1  CWS2«(CWATER/CSHEAR)**2 
SIG22«CWS2-1* 

XABwX* (CWS2/2#"1»“X**2)**2 

XB6*X* ( 1 • *X**2) *SQRT ( (CoSAL#*2*X**2) * ( SIG22»X**2) ) 
XC8«B*CWS2*#2*SQRT  (COSAL##2*X**2)  /A, 

C 

FA0.SQRT  <X*«2*C0SAL**2) *XBB/ < (XAB.XCB) **2.XBB**2) 

F8B«SQRT< (SORT! <X**2*P(lln**2*P(12) ).X**2-P<13) )/<  <X**2*P(11) )*•? 
1*P ( 12) ) ) 

C  P  ARRAY  CALCULATED  IN  SUBROUTINE  STPWP 
0NE1«FAB*F8B 
C 

99  RETURN 
END 


C*****FUNCTION  TWO**** 

C 

C 

C  CAGNIARD-ROSEnBAUM  INTEGRAnO  for  SLOW  NON-RIGID  BOTTOMS 

C 

C 

FUNCTION  TWO(X.P) 

DIMENSION  POO) 

COMMON  B.C0SALtC0STH,R2tSlN8E,SlNTH,CwATER»CB0T,CSHEAR»CST0N,RESID 
COMMON  C2.CB0T2.CSHR2,CBSH.C2SHR2,C*CR,SINTH2 
C 

SIGM2» (CWATER/CBOT ) **2-1 « 

FA0«X*SQRT  <SlGM2-X**2)/<  < 1,-B**2)*X**?*SIGM2«B**2) 

FBBaSQRT  < ( SORT <  <X**2*P < 1 1 ) j **2*p (12) ) ♦X**2-P ( 1 3) ) / (IX«*2«P ( 1 1 ) ) **2 
1*P ( 12) ) ) 

C  P  array  CALCULATED  IN  SUBROUTINE  STPWR 
TWO«FAB*FBB 

c 

RETURN 

end 


C ***** FUNCTION  SEVEN**** 

C 

C 

c  INTEGRAND  OF  CAGNIARD  INTEGRAL  USING  COMPLEX  ARITHMETIC. 

C 

C 

FUNCTION  SEVEN(Z.P) 

DIMENSION  P(30) 

COMMON  BiCOSAL iCOSTH»R2,SINBE,SINTH»CWATER»CBOT*CSHEAR»CSTON»RESID 
COMMON  C2.CBOT2.CSHR2.CHSH.C2SHR2.C4CB.SINTH2 
COMPLEX  F.RCOE.Y1.Y3.V.RT1.RT2.RT5.W.U1.U2.U3.U4 
COMPLEX  V2.XY1.XW 


ONE1001 
ONE1002 
ONE1003 
ONE1004 
ONElOOS 
ONE1006 
ONE1007 
ONE  10 08 
ONE  10 09 
ONE  1 0  1 0 
ONE  1011 
ONE1012 
ONE1013 
ONE  1014 
ONE  1015 
ONE1016 
ONE1017 
ONE101B 
ONElO 19 
ONE1020 
ONE1021 
ONE  1022 
ONE1023 
ONE1024 
ONE1025 
ONE  1026 


TWO  001 
TWO  002 
TWO  003 
TWO  004 
TWO  005 
TWO  006 
TWO  007 
TWO  O&B 
TWO  009 
TWO  OlO 
TWO  OH 
TWO  012 
TWO  0 1 3 
TWO  014 
TWO  015 
TWO  016 

two  oir 

TWO  018 
TWO  019 
TWO  020 


SEVN001 

SEVN002 

SEVN003 

SEVN004 

SEVN005 

SEVN006 

SEVN007 

SEVN008 

SEVN009 

SEVNOlO 

SEVN011 

SEVN012 
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C  P  ARRAY  CALCULATED  IN  SUBROUTINES  STPwA  AND  STPWB 
V-ChPLX  <P  < 1 ) • 7 ) 

V2«V*V 

RT1-CSGRT ( 1 .♦ v2) 

RT2«CSGRT(V2*CB0T2*C2) 

IF (CSHEAR.GT.O.)  GO  TO  20 
C 

c  non-rigid  bottoms 

Y3«0/C0OT*RT2 

RC0E«<PT1-Y3)/<RT1*Y3) 

GO  TO  40 

V- 

C  RIGID  BOTTOMS 

20  XY1«C2SHR2#V2*C2 

Y1«HT1«(XY1#XYUCBSH#V2#RT2*CSQ«T(V2*CSHR2*C2) ) 

Y3«C4C8*RT2 

RCOF  *<Yl-Y3)/(Yl*Y3) 

40  IF (7,EG,P  (A) )  GO  TO  50 

C  integrand  calculation 

XW«P ( 2) »C0STH*RT1 
X*V2#SINTH2*XW»X»( 

F-V/RTl/CSORT 
SEVEN-REAL (F* (RC0E-RT5) ) 

RETURN 

C 

C  CALCULATION  OF  ANS2 

50  RT5-RC0E 

U1«CMPLX<P(1) ,P(3) ) 

U2»1«*U1*U1 
U3-CS0RT (U2) 

XBa-P (2)*COSTH 

U4aHT5«CL0G( <  RT 1  *XB)  /  <CS<JRT  <U2«2. #XR*lj3*P  { 2)  **2-  SINTH2  WU3*XB) ) 
SEVEN«AIMAG<U4) 

RETURN 

End 


SEVN0I3 

SEVNOU 

SEVN015 

SEVNOU 

SEVN0I7 

SEVN0I8 

SEVN0I9 

SEVN020 

SEVN021 

SEVN022 

SEVN023 

SEVN024 

SEVN025 

SEVN026 

SEVN027 

SEVN028 

SEVN029 

SEVN030 

SEVN031 

SEVN032 

SEVN033 

SEVN03A 

SEVN035 

SEVN036 

SEVN037 

SEVN038 

SEVN039 

SEVN040 

SEVN041 

SEVN042 

SEVN043 

SEVN044 

SEVN045 

SEVN046 

SEVN047 


(^••••SUBROUTINE  EXE1****  EXE  1001 

C  EXE1002 

C  EXE1. 003 

C  CALCULATION  OF  EXPONENTIAL  INTEGRAL  El  TIMES  EXP  C-Y)  ,  NEGATIVE  Y  EXE1004 
C  EXE1005 

SUBROUTINE  EXEl(Y.ANS)  EXE1006 

DIMENSION  A(4)»B(4)*C<6)  EXE1007 

DATA  A/  8,5733287 ,18,059 017,8, 6347609,0,26777 373  /  EXE10C8 

DATA  B/  9, 5733223, 25. 632956,21, 0996531.3, 9584969  /  EXE1009 

DATA  C/  -0.57721566,0,99999193,-0,24991055,0.05519968,-0.00976004  EXE1010 

1  ,0.00107857  /  EXE1011 

X«-Y  EXE1012 

IF (X, IT, 1,0)  GO  TO  10  EXE1013 

ANS>(A(4)+X*(A(3)*X*(A(2)+X*(A(1)*X))))/(B(4)*X*(B(3)*X*(B(2)+  EXE1014 

1  X* (B  ( l ) 4X ) ) I ) /X  EXE1015 

RETURN  EXE1016 

10  ANS=F.XP(X)*(C(1)4X*(C<2)4X*(C(3)4X*<C(4)4X*(C(5J  *X*C(6)  )  ))  )-  EXE1017 

1  ALOG(X))  EXE1018 

RETURN  EXE1019 

END  EXE1020 
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c##«##subroutine  EXFI#### 

c 

c 

C  CALCULATION  OF  EXPONENTIAL  INTEGRAL  El  TIMES  EXP(-Y)  ,  POSITIVE  Y 
C 

SUBROUTINE  EXEI (Y.ANS) 

DIMENSION  P ( 1 0 ) *  A ( 6 ) 

EXTERNAL  EXPO 

DATA  A/  .25, .05555556, .01041667, .00166667, .00023148,. 00002834  / 

IF (Y.GT .0.5)  GO  TO  10 

U=Y*(1,+Y*(A(1)*Y*(A(2)+Y#(A(3)+Y*(A(A)*Y*(A(5)+Y*A(6)) )) > ) ) 
ANS«(0.57721566*ALOG(Y)+U)*EXP(-Y) 

RETURN 
10  P(1)=Y 

ANS1*FGI <l.,Y,4,EXPO,P) 

ANS«ANS1+1. 8951178  *EXP(-Y) 

RETURN 

END 


EXE  1 00 1 
EXE  1002 
EXE 1003 
EXEI00* 
EXEI005 
EXEI006 
EXE 1 007 
EXEI008 
exeiooo 
EXEIOlO 
EXE I  Oil 
EXE 1012 
EXEI013 
EXEIOU 
EXE 1 0 15 
EXE 1016 
EXEIOU 
EXE 1 0  18 


c*##««functICn  expo**** 

c 

c 

c  integrand  of  exponential  integral 
c 

FUNCTION  EXPO(X,P) 

dimension  p < i o i 
EXPO«EXP(X-P(l) )/X 

C  P(l)  IS  THE  ARGUMENT  OF  THE  EXPONENTIAL  INTEGRAL 
RETURN 
END 


EXPOOOl 

EXP0002 

EXP0003 

EXP0004 

EXPOOOS 

EXP0006 

EXP0007 

EXP0008 

EXP0009 

EXP0010 

EXPOOU 


C***«*FUNCTION  FGI**** 

C 

C 

C  THIS  SUBPROGRAM  INTEGRATES  THE  FUNCTION  F  BETWEEN  THE  LIMITS 

C  A  AND  B  USING  A  FOUR-POINT  GAUSSIAN  QUADRATURE  IN  EACH  OF  THE 

C  K  SUpINTERVALS.  P  IS  AN  ARRAY  USED  TO  TRANSFER  PARAMETERS  TO  THE 
C  FUNCTION  F, 

C 

FUNCTION  FGI (A,B,K,F,P) 

DIMENSION  V (4) »W(2) »SUM(4J ,P(1) 

DATA  V/  -.B61136311594053,-, 339981043584856, 

1  .3399810*358*856, .861136311594053  / 

DATA  W/  .347854845137454, .65214515*662546  / 

sum  < i >  «o • o 

SUM(2)=0.0 
sum (3) «o. o 
SUM (4) =0 .0 
H*(B-A)/FLOAT(K) 

H2*>H/2. 

AA«A+H2 
DO  20  L=1,K 
DO  10  Ini,* 

X«H2*V ( I )  4 AA 
10  SUM(I)cSUM(I)*F(X,P) 

20  AA=AA,H 

FGI*H2*(W(1)«(SUM(1)*SUH<*) ) +W (2) * (SUM (2) *SUM (3) ) J 

RETURN 

END 


FGI  001 
FGI  002 
FGI  003 
FGI  00* 
FGI  005 
FGI  006 
FGI  007 
FGI  008 
FGI  009 
FGI  010 
FGI  Oil 
FGI  012 
FGI  013 
FGI  014 
FGI  015 
FGI  016 
FGI  017 
FGI  018 
FGI  019 
FGI  020 
FGI  021 
FGI  022 
FGI  023 
FGI  024 
FGI  025 
FGI  026 
FGI  027 
FGI  02'3 
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SUBROUTINE  ploti***# 


plotting  subroutine  which  generates  a  plot  tape  for 
calcomp  plotter 

SUBROUTINE  PL0T1  <XX,YY, IPMAX.Xl *X2.ADATE*THE,WCH,CB0T,P0ISR*Z1) 
0IMENST0N  XX (1000) .YY(lOOO) ,BC0X<2> «BCDY(2) , TITLED (2) ,tITLE1 (2)  , 
1TITLE2 (2) t TITLE* (3) .TITLES (2) ,TITLE6 (2) .TITLE7 (2) ,TlTLEB(2) , 
2TITLE9<2) 

DATA  RCOX/10HTIME  (MICRt 10HOSEC)  /♦ 

1BCDY/10HPRESSURE  (*10HPSI)  /, 

2TITLE1 /lOHBOTTOH  REF.IOHLECTION  /, 
3TITLE2/10HCAGNIARD-R*10HOSENBAUM  /» 

*TITLE3/10HPLANE  WAVE/* 

STITLE4/10HPLANE  WAVE*10H  USING  CON,10HV.  INT.  /» 
6TITLE5/10HCOMPLEX  AR,10HITHHETIC  /» 

7TITLE6/10HANGLE  OF  I  * 10HNCIDENCE  /* 

BTITLE7/10HCHARGE  WEI.10HGHT  (LBS)  /. 

9TITLE8/1 OHCBOT  <FT/S»10HEC)  / 

DATA  TITLF9/10HPOISSON  RAilOHTIO  /» 

1TITLEO(1)/10HOATE  / 

TITLE0(2)*ADATE 
YMAX«#>.*X1 
YMlNa-3.#Xl 
YLIN«(YMAX-YMIN)/X1 
XLMAX«90 • 

XLHINsl . 

IYYMAX«0 

DO  *  1=1* I PMAX 

IF(YYd)-YMIN)  1*1*2 

1  YY(I)*VMIN 
GO  TO  * 

2  IF ( YY ( I ) “YMAX) 4*3*3 

3  YY(I)«yMAX 
I YYMAXal 

4  CONTINUE 
XMIN*XX  (1) 

XMAXsXX ( I PMAX ) 

CALL  SCAL  (X2.XMIN.XMAX,XLMIN,XLMAX»XLIN) 

IF (XMIN)5*10*10 

5  IF (XMAX) 10.10.6 

6  XS=XMIN 
XN»I. 

7  IF(ABSIXS)  l.E-36)9,9,71 
71  IF (XS)R*9*100 

a  xs«xmin*x2*xn 

XN*XN*1, 

GO  TO  7 
9  XN»XN-1,5 
YN«5.16 

CALL  CALCM1 (IPMAX.XX.YY.0.XMIN.XMAX*YMIN,YMAX.XLIN.YLIN.TITLE»0.BC 
1DX.15.BCDY.0 .FLOAT *18) 

CALL  SYMBL4(XN.YN,.14,BCDY,90,*1*) 

GO  TO  11 

10  CALL  CALCM1 (IPMAX.XX.YY. O.XMIN. XMAX. YMIN.YMAX.XLIN.YLIN. TITLE. O.BC 
idx.i5.bcoy.i4, float, is) 

11  IF(IYYMAX) 100,12*13 

12  XS»(XMAX-XMIN)/X2-3. 

GO  TO  14 

13  XS«(XX (IYYMAX)-XMIN)/X2*1. 

14  CALL  SYMBL4  (XS, 9*. .14* TITLE  1*0, *20) 

15  IF (Z1 -1 • ) 1 7 , 18* 16 

16  IF(Zl-3.) 19,20*100 


PLOTOOl 
PLOT002 
PLOT003 
PLOT004 
PLOT005 
PLOT006 
PLOT007 
PLOT008 
PLOT009 
PLOTOIO 
PLOTOU 
PLOT012 
PLOT013 
PLOT014 
PLOT015 
PLOT016 
PLOT017 
PLOT018 
PLOT019 
PLOT020 
PLOT021 
PLOT022 
PLOT023 
PLOT024 
PLOT025 
PLOT026 
PLOT027 
PLOT028 
PLOT029 
PLOT030 
PLOT031 
PLOT032 
PL0T033 
PLOT034 
PLOT035 
PLOT036 
PLOT037 
PLOT038 
PLOT039 
PLOT040 
PLOT041 
PLOT042 
PLOT043 
PLOT044 
PLOT045 
PLOT046 
PLOT047 
PLOTO40 
PL0T049 
PLOT050 
PLOT051 
PLOT052 
PLOT053 
PLOT054 
PLOT055 
PLOT056 
PLOT057 
PLOT058 
PLOT059 
PLOT  060 
PLOT061 
PLOT062 
PLOT063 
PLOT064 
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17  CALL  SYMBL4  <XS,8.?,.14.TITLE2.0.,20) 

18  «LL°sviBL4(XS. 8.7. .14, TITLES, 0. .10) 

19  CALL0SYMBL4  (XS,8.7,.14,TITLE4,0.,30) 

„  °°J°CvisL4  (XS,8.7,. 14, TITLES, 0. ,20) 

?!  "tt  SVMBL*  (XS,B.*..1*.TITLE6.0..?»I 

XN*XS+. 14*18. *.2  rue . n  .2) 

CALL  NUMBR(XN,8.4,.14,THE.0.,2 

call  SYMBL4  (XS»8,1 , .14 ,TITLE7,0. 

XN*XS*.14*19.*,2  5> 

CALL  NUMBR (XN,8. 1 , .14 ,WCH»0* ,5) 

CALL  SYMBL*  ( XS,7 .8, . 14, T ITLES, 0 . ,2« ) 

XN»XS*.14*13,*.2  ranT.n. .2) 

S!± 

CALL  SYMBL4  (XS, 7. 2, .14, TITLED, 0. ,2 
CALL  CALCMKO.O.) 

RETURN 

100  WRITE (6,22) 

CALL  CALCMKO.O.) 

2a  FORMAT 11H1 »10X«3*HPL0,TING  ERROR  IN  SOBROUT.BE  RUOTl  //« 

END 


PLOT  065 

PLOT066 

PLOT067 

PLOT  068 

PLOT069 

PLOT070 

PL0T071 

PL0T072 

PLOT073 

PLOT074 

PLOT075 

PLOT076 

PLOT077 

PLOT078 

PLOT079 

PLOTOSO 

PLOT081 

PL0T082 

PLOT083 

PLOT084 
PLOT085 
PLOT  086 
PLOT087 

PLOT088 

PLOT089 

PLOT090 

PLOT091 


C««****SUBR01JTINE  SCAL**** 

SUBROUTINE  FOR  SCALING  PLOTS 

SUBROUTINE  SCAL  (XSCALE.XMIN,XMAX,XLMIN,XLMAX,XLIN) 

DIMENSION  X (6) 

SXMAX=XMAX 

SXMIN*XMlN 

1  IJ*1 

2  Ir  (XSCALE)3»3«6 

c  determination  OF  SCALE  RANGE 

C  3  M1*ALOG10 (ABS (XMIN) ) 

4  XMsMl 
XMslO.**XM 
GO  TO  7 

5  XM*M1-1 
XM=10.**XM 
GO  TO  7 

6  XMsXSCALE 

C  SCALE  FACTORS  ALLOWED 

°  7  X(1)*1.*XM 

X(2)=2.*XM 
X (3) *2 ,5*XM 
X (4 ) »5  »*XM 
X (5) *7 ,5*XM 
X (6) si ft.*XM 


SCAL001 

SCAL002 

SCAL003 

SCAL004 

SCAL005 

SCAL006 

SCAL007 

SCAL008 

SCAL009 

SCALOlO 

SCAL011 

SCAL012 

SCAL013 

SCAL014 

SCAL015 

SCAL016 

SCAL017 

SCAL018 

SCAL019 

SCAL020 

SCAL021 

SCA1.022 

SCAL023 

SCA1.024 

SCAL025 

SCAL026 

SCAL027 

SCAL028 

SCAL029 

SCAL030 

SCAL031 

SCAL032 

SCAL033 
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IF(XSCALE)8,B,24 
AUTOMATIC  SCALING 


8  DO  21  I  »1 ,6 

DETERMINATION  OF  SCALED  MINIMUM 

XMINsSXMlN 

XMAX=SXMAX 

NSCALE=XMIN/X<I) 

IF (NSCALE >11,9,12 

9  IF(XMIN) 11.13,10 

10  XMINaO. 

GO  TO  13 

11  XNsNSCALE-1 
XMIN=XN*X ( I ) 

GO  TO  13 

12  XN=NSCALE 
XMIN-XN*X(U 

DETERMINATION  OF  SCALED  MAXIMUM 

13  NSCALE»XMAX/X(U *1. 

XNsNSCALE 

XMAX»XN*X(I) 

LENTH  OF  SCALE  AXIS  CALCULATED  AT  THIS  POINT 

XLIN=(XMAX-XMIN)/X(I) 

IF(XLIN-XLMAX) 14.18.17 

14  IF(XLlN-XLMIN)15.18,ia 

15  IF  (1-1)35,16.16 

16  XNsXM*l.E*01 
IJ*IJ+I 

IF ( I J-4) 7 .7,23 

17  IF  (1-6)21,177,35 
177  XMaXMM.E-01 

I J“I J+l 

IF ( I J-4) 7,7,23 

18  IF (XSCALE) 19,20,35 

19  IF(ABS(XSCALE)-X(I))22,20,22 

20  XSCALEsX(I) 

GO  TO  39 

21  CONTINUE 

22  XSCALEsABS (XSCALE) 

XMlNaSXMlN 

XMAXsSXMAX 

GO  TO  1 

23  IF (XSCALE)38.36,35 
FIXED  SCALIN' 

24  DO  25  I»1 ,6 

IF (XSCALE-X(I) >25.26,25 

25  CONTINUE 
GO  TO  37 

DETERMINATION  OF  SCALE  MINIMUM 

26  NSCALE=XMIN/XSCALE 
IF (NSC ALE) 28,27,30 

27  IF ( XMIN) 28 ,31 ,29 


SCAL034 

SCAL035 

SCAL036 

SCAL037 

SCAL038 

SCAL039 

scaloao 

SCALC41 

SCAL042 

SCALC43 

SCAL044 

SCAL045 

SCAL046 

SCAL047 

SCAL048 

SCAL0'9 

SCAL050 

SCAL0S1 

SCAL052 

SCAL053 

SCAL054 

SCAL055 

SCAL056 

SCAL057 

SCAL058 

SCAL059 

SCAL060 

SCAL061 

SCAL062 

SCAL063 

SCAL064 

SCAL065 

SCAL066 

SCAL067 

SCAL068 

SCAL069 

SCAL070 

SCAL071 

SCAL072 

SCAL073 

SCALC74 

SCAL075 

SCAL076 

SCAL077 

SCAL078 

SCAL079 

SCAL080 

SCAL081 

SCAL082 

SCAL083 

SCA1.084 

SCAL085 

SCAL086 

SCAL087 

SCAL088 

SCAL089 

SCAL090 

SCAL091 

SCAL092 

SCAL093 

SCAL094 

SCAL095 

SCAL096 
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28  XN=NSCALE-1  SCAL097 

XMlN»XN*XSCALE  SCAL098 

GO  TO  31  SCAL099 

29  XMlNsO.  SCAL100 

GO  TO  31  SCALlOl 

30  XN*NSCA1_E  SCAL102 

XMINsXN*XSCALE  SCAL103 

SCAL104 

DETERMINATION  OF  SCALED  MAXIMUM  SCAL105 

SCAL106 

31  NSCALE=XMAX/XSCALE*1.  SCAL107 

XNsNSCALE  SCAL108 

XMAX=XN*XSCALE  SCAL109 

SCAL110 

LENGTH  of  scale  AXIS  SCAL111 

SCAL112 

XLlNa ( XMAX-XMIN; /XSCALE  SCAL113 

IF(XLIN-XLMAX)32*39*33  SCAL114 

32  IF (XLlN-XLMIN) 34*39*39  SCAL1 15 

33  XSCALEaXSCALE«l.E*01  SCAL116 

IJ=IJ*1  SCAL117 

IF(IJ-4)2,2,37  SCAL118 

34  XSCALE*XSCALE#1.E-01  3CAL119 

IJCIJ*1  SCAL120 

IF(IJ-4)2,2,37  SCAL121 

35  STOP  SCAL122 

36  WRITE(6,200)XMAX.XMIN.XLIN, <X(I)»Isl»6>  SCAL123 

GO  TO  35  SCAL124 

37  WRITE(6,201JXMAX.XMIN,XLIN,  (X(I)  »Iel,6)  SCAL125 

GO  TO  35  SCAL126 

38  WRITE(6t202)XMAX.XMIN,XLlN, (X  < I ) #I»1»6>  SCAL127 

GO  TO  22  SCAL12B 

200  FORMAT (42X*36HAUT0MATIC  SCALING  CANNOT  BE  ACHIEVED/30X*5HXMAX*1E14  SCAL129 

1.5.2X,5HXMIN*lEl4.5.2X.5HXLINclE14,5//51Xfl7HSCALE  FACTORS  ARE//18  SCAL130 
2X.6E14.5//)  SCAI.131 

201  FORMAT (45X«32HFIXED  SCALING  CANNOT  be  ACHIEVED/30X*5HXMAX*lE14.5t2  SCAL132 

1X,5HXMIN=1E14.5,2X.5HXLIN=1E14.5//51X,17HSCALE  FACTORS  ARE//18X,6E  SCAL133 
214,5//)  SCAL134 

202  FORMAT (9X.76HAUT0MATIC  SCALING  CANNOT  ACHIEVE  DESIRED  SCALE  FACTOR  SCAL135 

1.WILL  TRY  FIXED  SCALING//30X»5HXMAX=1E14,5.2X,5HXMIN=1E14.5.2X»5HX  SCAL136 
2LIN=lE14»5//5l Xt 17HSCALE  FACTORS  ARE//18Xf 6E14.5//)  SCAL137 

39  RETt'HN  SCAL138 

ENl,  SCAL139 
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C  #####  pTv  PROGRAM  ••••• 

c 

SUBROUTINE  PTV<TImER2,T3,T4,T5.RAD»PTS*0PTI0N.C0SA,RH0W»CWAT» 

1  T.V.VS) 

C 

U  THIS  subprogram  controls  the  iteration  for  the  peak 
c  translational  velocity*  PTV.  it  is  the  only  subroutine  of  the 

C  PTV  PROGRAM  WHICH  IS  CALLED  FROM  ThE  MAIN  PROGRAM. 

C 

DIMENSION  QX ( 1 000 ) tQY(lOOO) *  IS (2) 

DIMENSION  6 (6) 

DIMENSION  A  (50 ) *C (50 ) 

COMMON  /GXY/QX.QY 
COMMON  /GIS/IS 
C 
C 

IF (OPTION.GT.O. )  GO  TO  10 
WRITE  (6*580) 

WRITE (6.600)  TIMER2. T3.T4.T5.RAD.PTS. OPTION. COSA.RHOW.CWAT 
WRITE (6.590) 

C 

C  74.21457  IS  A  UNITS  CONVERSION  FACTOR 
10  VC=2.#74.21457/RH0W/RA0 

n=pts 

T=T4 

DT=(T5-T)/Fl0AT(N-1) 

IF(T.LE.O.)  N=N-1 
IF(T.LE.O.)  T*DT/2. 

IS  < 1 ) *2 

IS'2)=1 

G(2)=CwAT/RAD 

G(3)=TIMER2 

G(5)=SQRT(TIMER2-T3) 

G(6)=SQRT(TIMER2) 

C  INITIAL  SEARCH  FOR  MAXIMUM  VELOCITY 

DO  40  1  =  1  *N 
G(1)xT 
V=VC*FV(G) 

A(I)=T 

C(I)=V 

VS=2.#C0SA«V 

IF (OPTION.LE.O. )  WRITE (6.610)  T.V.VS 
T=T*DT 
40  CONTINUE 

C  ITERATION  FOR  PTV 

C  DETERMINE  THE  MAXIMUM  VELOCITY  FROM  C  ARRAY 
CALL  XMAX (C.NyM.Ml ) 

A2=A  (Ml ) 

C2=C (Ml ) 

A(1)*A(M) 

C ( 1 ) =C (M) 

A (2) =A2 
C (2) =C2 
DA=DT 

T=A ( 1 ) — 1 • 8  *D  A 
IF(T.LE.O.)  T*DA/5. 

UT=DA/2. 

DO  45  1=3.10 

G(1>=T 

V=VC*FV(G) 

A ( I ) =T 
C(I)=V 

VS=2.*C0SA*V 

IF(OPTION.LE.O.)  WRITE<6.610)  T.V.VS 
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T=T  *DT 
45  CONTINUE 
Nslo 

IE (IABS (M-Ml ) .LT.3)  GO  TO  55 

T=A(2)-0.8*DA 

IEU.LE.O.)  TrDA/5. 

0T»0A/3. 

DO  50  1=11.16 

G(1)=T 

V=VC*FV(G) 

A(I)=T 

C(I)=v 

VS=2.*C0SA*V 

IF(OPTION.LE.O.)  WRITE (6.610 )  T.V.VS 
T=T +OT 
50  CONTINUE 
N=l6 

55  CONTINUE 
DO  75  JJ=1 .6 
CALL  XMAX(C.N.M.Ml) 

IE ( JJ.LT #3)  GO  TO  62 

IE  (ABS ( (C(M)-C(Ml) )  /C  (M) ) *LT. 0*001)  GO  TO  110 
IE (JJ.EQ.6)  GO  TO  120 
62  Nalo 
T1=A(M) 

T2=A (Ml ) 

VI =C (M) 

V2=C (Ml ) 

A(9)=T1 
A  ( 10 ) =T2 
C(9)=V1 
C ( 1 0 ) =V2 

DT=ABS(Tl-T2)/5. 

11  =  1 

DO  70  1=1,8 

T=Tl+OT*FLOAT( (I-10)/2«II) 

IE  (T  ,L£.  0.0)  GO  TO  64 

G ( 1 ) =T 

V=VC*FV(G) 

VS=2,*COSA*V 
GO  TO  66 

C  WHEN  T  IS  LESS  THAN  ZERO  SET  TO  ZERO, 

64  T  =  0.0 
V  =  0.0 

66  *F  (-PTI°N  »<-E.  °.°>  “RITE  (6,610)  T.V.VS 

C(I)=V 
II=-1*II 
70  CONTINUE 
75  CONTINUE 
110  V=C(M) 

T=A (M) 

VS=2.*C0SA*C(M) 

IE (OPTION, LE.O. )  WRITE(6.620)  A(M),C(M)»VS 
RETURN 
120  V=C (M) 

T=A(M) 

VS=2.*C0SA*C(M) 

VSl='j.*COSA*C(Ml) 

WRITE (6.630)  T. V.VS.A (Ml ) ,C (Ml ) , VS1 
RETURN  1 


580  FORMAT (1H1.1 OX. 30HTRANSLATIONAL  VELOCITY  PROGRAM  ) 
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590  FCRMAT(1HO*5X,45HITERATION  FOR  PEAK  TRANSLATIONAL  VELOCITY  PTV  // 

1  12X»9HTlME (SEC) »8X» 16H VELOCITY (FT/SEC)  .3X.25H VERTICAL  VELOCITYIF 
2T/SEC)  /29X.16HTARGET  SUBMERGED. 7X,17hTARQET  AT  SURFACE  ) 

600  FORMAT ( 1 HO  »5X ,23H INPUT  TO  SUBROUTINE  PTV  //  10X. 

1  45HTIMER2.T3.T4.T5, RAD. PTS. OPTION, COSA.RHOW.CVAT  //1P5E14.5/ 

2  1P5E14.5  ) 

610  FORMAT (1P3E22, 6) 

620  FORMAT(1HO,6X,20H*********»»»»**»»«»*  .9X.3HPTV, 19X.3HPTV// 

1  1P3E22.6) 

630  FORMAT (1H0,42H**»  WARNING  ITERATION  DID  NOT  CONVERGE  •**  ,5X. 

1  35n"AXlMUM  ANO  NEAREST  VALUE  ARE  RIVEN  // 

1  12X.9HT1ME (SEC) »8X, 16HVEL0CITY (FT /SEC)  .3X.25HVERTICAL  VELOCITYIF 
2T/SEC!  /29X.16HTARGET  SUBMERGED, 7X, 17HTARGET  AT  SURFACE  / 

3  (1P3E22.6)) 

END 


FUNCTION  FV (G) 

THIS  SUBPROGRAM  SETS  UP  THE  INTEGRATION  FOR 
THE  TRANSLATIONAL  VELOCITY  V 

DIMENSION  G (6 ) 

EXTERNAL  FI 
DATA  N/18/ 

NN=FLOAT(N)*G(l)*G(2)/B. 

NNrMAXO (NN.fl) 

NNsMINO (NN.N) 

X*G ( 1 ) -8./G (2) 

IF  (X.6T«G(3) )  GO  TO  43 
Z1*G<6> 

IFCX.GT.O.)  Z1»SQRT (G(3)-X) 

IF (G ( 1 ) »GT »G (3 ) )  GO  TO  40 
G (4) s-1 , 0 

Z2*SQRT (G (3) -G ( 1 ) ) 

INTEGRATION  FOR  T  .LE.  TIMER2 
FVr-FGI (ZI.Z2.NN.F1.3) 

RETURN 
40  Z2=0. 

Z3»SQRT<G(1)-G(3) > 

IF(G(3).EQ.O.)  GO  TO  45 
G<4)«-1.0 

N1«Z1/(Z1*Z3) *FLOAT  (NN) ,2.0 
NNN«Z3/(Z1.Z3)»FLOAT<NN).2,0 

INTEGRATION  FOR  INTERVAL  WHICH  INCLUDES  TIMER2 
Vls-FGKZ1.Z2.N1, FI. G) 

GUjsl.O 

V?=FGI (Z2.Z3.NNN. FI, G) 

FVSV1+V2 

RETURN 

43  Z2»SQRT(X-G(3) ) 

Z3»SQRT (G ( 1 ) «G (3) ) 

45  G<4)*1.0 

INTEGRATION  FOR  T  LARGER  THAN  TIMER2  BUT  THE 
INTERVAL  DOES  NOT  INCLUDE  TIMER2. 

FVsFGI (Z2,Z3,NN,Fl«G) 

RETURN 

END 
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FUNCTION  Fl(Z.G) 

THIS  SUBPROGRAM  CALCULATES  THE  PRODUCT  INCIDENT  PRESSURE  * 
REDUCED  STEP  WAVE  ACCELERATION  BY  CALLING  THE  INTERPOLATION 
PROGRAMS  VTAB  AND  PTAB. 

DIMENSION  QX(IOOO) »QY(1000) .IS (2) 

DIMENSION  G{6) *QQX ( 120) »QQY ( 120) 

COMMON  /QXY/QXtQY 
COMMON  /QIS/IS 

REDUCED  STEP  WAVE  ACCELERATION  OF  a  CYLINDER 

DATA  (QQX(I) *1=1*106)  /0. * . 0125, .025, .0375 , .050 . . 075. . 100 . 

1  .125., 150* .175,. 200*. 225,. 250*. 275*. 300,. 325,. 350*. 375, 

2  .4000,, *25,. *50.. 475*. 500,. 525,. 550,. 575*. 600*. 625. .650. 

3  .675*. 700,. 725,. 750,. 775,. 800,. 825,. 850,. 875,. 900,. 925* .950* 

4  .975,1.00,1.05,1,10,1.15,1.20,1.25,1.30,1.35,1.40,1.45, 

5  1.50, 1.55, 1.60, 1.65, 1.70, 1.75, 1.80, 1.85, 1,90, 1.95, 2. On. 

6  2, 05, 2, 10, 2, 15, 2‘.20, 2. 25, 2. 30, 2. 35, 2. 40, 2. 45, 2. 50, 2, 55, 

7  2. 60, 2. 65, 2. 70, 2. 75, 2. 80, 2.85:2.90,3. 00, 3, 10. 3.2, 3. 3, 3. 4, 

8  3. 5, 3.6, 3. 7, 3. 8, 3. 9, 4.0, 4, 2, 4. 4, 4, 6, 4. 8, 5.0, 5. 25, 5, 50. 

9  5. 75 .6.00. 6. 25 ,6.5. 7, 0.7. 5, 8. j  / 

DATA  (QQY(I) .1=1.60)  /  0*0,  .198193.. 275935,. 332694., 378180, 

1  . 448836,. 502189,. 544000*. 577342,. 6041 11,. 625589,. 642701, 

2  .656143,. 666457,. 674079 ,.679365,. 682612 ,,684070,. 683955, 

3  .682452, .679721 , .675904, .671127, .665499, .659120, .652078, 

4  .644453, .636315, .627730, .618755, .609444 , .599844, .589999, 

5  .579949, .569730, .559374, .548913, ,538372, .527777, ,517151 , 

6  .506515, .495887, .485284, .464215, .443417, .422977, ,402968, 

7  .383447 , ,364460 , .3460*2, .328218, .311 008 , .294424, ,278471 , 

8  .263152, .248465, .234404, .220960, .208124,. 195881  / 

DATA  (QOY(I), 1=61, 106)  /  . 184219, . 173122, .162573, . 152555, 

1  . 143051,. 134041,, 125509,. 1 17435, -109801,. 102590,. 095782, 

2  .089361,. 083308,. 077608,. 072242,. 067196 , .062453, , 057999, 

3  .053818,. 049897,. 046221,. 039556,. 033725,. 028637,. 024209, 

4  .020368,. 017044,. 014177,. 01 17 12,. 009599,. 007795,. 006260, 

5  .003863,. 002172,. 001009,. 000230, -0.000267,. 0.000619, 

6  -0.000774,-0.000804,-0.000767,-0.000696,-0,000606, 

7  -0.000430,-0.000297,-0.000206  / 

IF (G(4).GT.O.)  GO  TO  20 
X*G (3) -Z*Z 
GO  TO  30 
20  X«G (3) +2*Z 
30  XD« (G ( 1 ) -X) *G (2) 

IF(Z.GT.G(5) )  GO  TO  35 
P=PTAB,'X*0X.0Y,IS(2)  ) 

00  TO  40 

35  P*VTAB (X,QX«QY«IS(2) ) 

40  Fl=Z#P#VTAB(XD,OOX,QQY,IS(l) ) 

RETURN 

END 
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SUBROUTINE  XMAX {B.N.M.Ml ) 

C 

c 

c  THIS  subprogram  determines  the  locations  of  the  two  largest 

C  ABSOLUTE  VALUES  OF  MEMBERS  OF  THE  B  ARRAY. 

C 

C 

DIMENSION  B (50 ) 

X»ABS (B ( 1 ) ) 

M=1 

DO  10  I>2«N 

IFtABStB(I) ) .LT.X)  GO  TO  10 

M»I 

X=ABS(B(M) ) 

10  CONTINUE 
Ml  =  1 

IF(M.EO.l)  Ml*2 
XsABS (B (Ml ) ) 

DO  20  I«2.N 

IF ( ABS (B ( I ) ) .LT.X)  GO  TO  20 

IF ( I iEO.M)  GO  TO  20 

M1»I 

X-«ABS(B(M1) ) 

20  CONTINUE 
RETURN 
END 


FUNCTION  VTAB<X,Y,Z.K) 

this  SUBPROGRAM  PERFORMS  A  SECOND  ORDER  LAGRANGlAN  INTERPOLATION 

THE  INDEPENDENT  VARIABLE  IS  STOREO  IN  THE  Y  array  in  increasing 
ORJER.  THE  DEPENDENT  VARIABLE  IS  STORED  IN  THE  Z  ARRAY. 

X  IS  THE  POINT  AT  WHICH  THE  FUNCTION  IS  TO  BE  EVALUATED. 

K  IS  THE  NUMBER  OF  THE  ELEMENT  IN  THE  V  ARRAY  WHICH  IS  FIRST 
COMPARED  WITH  X. 


DIMENSION  Y (1000) *Z  (1000) 

IFIX.LE.O.)  GO  TO  50 
DO  10  I»K* 1000 
J=I 

IF(Y{I).GT.X)  GO  TO  20 
10  CONTINUE 
20  J*MAX0(3#J-1) 

DO  30  IaltlOOO 

IF (Y (J) .LT.X)  GO  TO  40 

jaj-l 

IFU.LT. 3)  GO  TO  40 
30  CONTINUE 
40  K»J* 1 

IF (Z ( J) .EQ.Z (K) )  GO  TO  60 
L*  J-I 

A»(X-Y(K))/(Y(J)-Y(L>) 

C«IX-YILI)/|V|K)-YIJ)I 

IF ( (A.LT.-5.0) .OR. (C.GT.5.0) )  GO  TO  60 

B*(X-Y(J))/(Y(K)-Y(L)) 

VTAB*C*(B*Z(K)-A*Z<J> )*A*B*Z(L> 

RETURN 
50  VTAB*0. 

RETURN 

60  VTAB«Z(J)*(X-Y(J))*(Z(K)-Z(J)  ) / ( V (K ) -Y ( J)  ) 
RETURN 
END 
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FUNCTION  PTAB  (XiY.ZtK) 


This  subprogram  performs  a  second  ORDER  lagrangian  INTERPOLATION 
WITH  PROVISIONS  FOR  HANDLING  A  SINGULARITY. 

FUNCTION  ARGUMENTS  ARE  THE  SAME  AS  IN  VTAB  . 


DIMENSION  Y  < 1000 ) .2  (1000) 

IF (X.LE.0. )  GO  TO  50 
DO  If  IsK.1000 
J=I 

IF(Yd).GT.X)  GO  TO  20 

10  continue 

20  J=MAX0 (3* J-l ) 

DO  30  IsltlOOO 
IF(Y(J).LT.X>  GO  TO  40 
J=J-1 

IF ( JtLT *3)  GO  TO  40 
30  CONTINUE 
40  J=J*1 
JJ*J 

THE  FOLLOWING  THREE  STATEMENTS  PROVIDE  FOR  EXTRAPOLATION 
AROUND  A  SINGULARITY. 

IF(ABS(Z(J) ) .GT.1.0E20) JJaJ-2 
IF (ABS  (Z (J-l) ) .GT.1.0E20) JJeJ+1 

IF  (  ( JJ.EQ.J)  .AND.  (ABS  (Z  (J-2) )  .LT.1.0E20)  )  JJ«=J-1 

jsjj 

K=J*1 

IF (Z( J) .EQ.Z (K) )  GO  TO  60 
L- J-l 

A=(X-Y (K) )/(Y (J)-Y(L* ) 

Cs(X-Y(L))/(Y(K)-Y«J)) 

IF ( (A.LT.-5.0) .OR. (C.GT.5.0) 1  GO  TO  60 
B*(X-Y(J) ) / ( Y  <K) -Y (L) ) 

PTAB«C*(B*Z(K)-A»Z(J) )*A«B#Z(L) 

RETURN 
50  PTABsO. 

RETURN 

60  PTAB*Z(J)*(X-Y(J) )*(Z(K)-Z(J) )/(Y(K)-Y(J> ) 

RETURN 

END 
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APPENDIX  B 

SAMPLE  PROGRAM  OUTPUTS 
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REDIJCEO  TIME  CONSTANT  OF  INCIDENT  WAVE  THETA-  .2709AA53E-0I 

ACTUAL  S“  TIME  CONSTANT  IN  MILLISECONDS  THET-  , 22783* 30E-0 I 

REDUCED  TIME  CONSTANT  OF  BOTTOM  REFLECTED  WAVE  THETAR-  .27I89556E-0I 

BOTTOM  PFFLECTED  WAVE  TIME  CONSTANT  IN  MILLISECONDS  THETR-  .22B6340IE-OI 
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B0TT0H  REFLECTION 
CAGNIAR0-R0SENBAUH 
ANGLE  OF  INCIDENCE  72.03 
CHARGE  WEIGHT  (LBS)  .00125 
CB0T  (FT/SEC)  5870.00 
POISSON  RATIO  .50000 
DATE  06/17/71 
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FIG .  2  SAMPLE  CALCOMP  PLOT  OF  PRESSURE-TIME  HISTORY 
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APPENDIX  C 

TWO  SPECIAL  INPUT  OPTIONS 

Two  special  input  options  for  altering  the  input  geometry  are 
provided  in  the  code  BOTREF  using  the  variables  THOVAL  and  Z2.  The 
first  option  allows  the  programmer  to  specify  an  incident  angle  6 
which  is  expressed  as  the  ratio  of  the  incident  angle  6  to  the 
critical  angle  9cr.  This  is  accomplished  by  setting  THOVAL  =  0/®cr>  °* 
If  THOVAL  *  0,  this  option  is  ignored.  The  programmer  also  must 
supply  the  water  depth  H,  the  horizontal  range  r,  and  the  gauge 
depth  dg.  If  possible,  the  code  calculates  the  required  charge 
depth  d  keeping  dg  fixed.  Otherwise,  d  is  set  to  H,  and  a  new 
value  of  dg  is  determined.  Thus  this  option  permits  the  user  to 
calculate  bottom  reflections  for  a  range  of  incident  angles  without 
first  having  to  determine  the  exact  geometry  that  is  required. 

The  second  option  using  Z2  >  0  provides  an  alternative  means 
of  specifying  the  reflection  geometry.  The  arrival  time  Z2  (in 
microseconds)  of  the  bottom  reflected  wave  after  the  direct  wave  is 
exceedingly  sensitive  to  the  geometry  which  often  cannot  be  measured 
with  the  necessary  accuracy.  This  time,  which  can  be  accurately 
measured,  provides  the  means  which  can  be  used  to  correct  the  input 
geometry.  The  geometry  is  changed  by  altering  d  and  dg  and  holding 
r  fixed  so  that  the  incident  angle  9  and  the  bottom  reflected  slant 
range  Rr  are  unchanged.  Typically,  the  direct  and  reflected  pulses 
are  only  slightly  altered,  but  the  change  in  their  sum  may  be  signifi¬ 
cant. 
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